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ABSTRACT 


This  thesis  develops  and  validates  a  satellite  trajectory  optimization  model.  A 
summary  is  given  of  the  general  mathematical  principles  of  dynamic  optimal  control  to 
minimize  fuel  consumed  or  transfer  time.  The  dynamic  equations  of  motion  for  a 
satellite  are  based  upon  equinoctial  orbital  elements  in  order  to  avoid  singularities  for 
circular  or  equatorial  orbits.  The  study  is  restricted  to  the  two-body  problem,  with  engine 
thrust  as  the  only  possible  perturbation.  The  optimal  control  problems  are  solved  using 
the  general  purpose  dynamic  optimization  software,  DIDO.  The  dynamical  model 
together  with  the  fuel  optimal  control  problem  is  validated  by  simulating  several  well 
known  orbit  transfers.  By  replicating  the  single  satellite  model,  this  thesis  shows  that  a 
multi-satellite  model  which  optimizes  all  vehicles  concurrently  can  be  easily  built.  The 
specific  scenario  under  study  involves  the  injection  of  multiple  satellites  from  a  common 
launch  vehicle;  however,  the  methods  and  model  are  applicable  to  spacecraft  formation 
problems  as  well. 
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I.  INTRODUCTION 


A.  PROBLEM 

In  the  space  business,  available  fuel  can  be  directly  equated  to  mission  life. 
Depletion  of  a  satellite’s  fuel  and  the  associated  loss  of  maneuvering  capability  generally 
renders  useless  whatever  remaining  capabilities  a  satellite  may  have.  Add  to  this  a 
tremendous  cost  of  launching  satellites,  where  the  mass  of  the  payload  needs  to  be 
maximized,  and  careful  trades  must  be  made  to  make  sure  every  ounce  of  fuel  onboard 
has  a  purpose.  Therefore  the  task  of  ensuring  efficient  maneuvering  during  the  life  of  the 
vehicle  takes  great  scrutiny,  for  every  gram  of  fuel  saved  early  is  a  gram  that  can  be  used 
later  to  extend  the  mission  life. 

Designing  fuel  efficient  maneuvers  for  a  single  vehicle  is  a  time  consuming 
process,  involving  iterative  checks  to  ensure  that  the  minimum  amount  of  fuel  is 
expended.  Attempting  to  do  this  concurrently  for  multiple  vehicles  is  very  difficult; 
instead,  a  more  serial  approach  from  vehicle  to  vehicle  is  usually  preferred.  For  some 
scenarios,  the  process  does  not  necessarily  take  into  account  relationships  and 
interactions  between  vehicles.  In  situations  where  spacecraft  are  closely  spaced,  parallel 
computations  are  difficult  to  perform  because  engineers  must  understand  where  the  first 
satellite  is  and  where  it  will  go  before  they  can  calculate  what  the  second  should  do. 
Therefore,  it  is  a  mathematically  complex  and  labor-intensive  process,  involving 
numerous  semi-automated  tools. 

These  methods  are  tried  and  true,  and  to  some  degree  are  considered  black  art 
practiced  by  highly  knowledgeable  and  experienced  specialists.  The  space  industry  is 
notorious  in  its  conservatism  and  caution  in  moving  away  from  “what  works”  for 
promises  of  improvement  that  don’t  always  bear  fruit,  but  one  can  still  question  whether 
these  methods  are  the  most  effective  use  of  time  and  energy  given  the  advancements  in 
computation  and  optimization.  An  automated  tool  that  could  handle  this  process  without 
fail  for  any  circumstance  would  be  ideal;  however,  this  is  a  highly  unrealistic  prospect. 
Perhaps  nearer  at  hand  are  methods  and  tools  that  can  quickly  provide  a  starting  point  to 
this  process,  enabling  a  much  shortened  timeline  to  completion.  It  could  also  enable 
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quick  calculation  of  several  different  possible  seenarios  to  choose  from,  with  a  minimal 
resource  impact.  The  research  for  this  thesis  explores  these  new  methods. 


B,  PROPOSED  SOLUTION 

Tools  have  been  developed  over  the  past  few  years  at  the  Naval  Postgraduate 
School  that  can  be  used  to  great  effect  in  almost  any  imaginable  dynamic  optimization 
problem.  Specifically,  a  generic  optimization  engine,  DIDO,  has  been  developed  whieh 
allows  the  user  to  quiekly  perform  dynamic  optimization  on  any  problem  that  can  be 
properly  framed  in  mathematical  terms,  and  this  tool  has  been  at  the  heart  of  several 
reeent  theses  and  dissertations  published  at  the  sehool  [King  (2002),  Stephens  (2002), 
Josselyn  (2003),  Shaffer  (2004),  Fleming  (2004)].  The  generic  capability  that  DIDO 
presents  is  startling,  for  the  only  resemblance  of  these  problems  (outside  of  the 
eoincidence  of  all  being  aerospaee  related)  is  in  the  underlying  DIDO/MATLAB^m 
toolset. 

The  primary  goal  of  this  researeh  is  to  show  that  a  DIDO  based  optimization  tool 
ean  be  developed  whieh  can  predict  fuel  optimal  control  maneuvers  for  orbital 
applieation.  Some  previous  work  has  been  done  in  this  area  speeific  to  single  satellite 
models.  However,  this  research  aims  to  extend  this  work  most  notably  by  extending 
developed  ideas  and  methods  to  optimal  control  prediction  for  multiple  satellites 
simultaneously.  Additionally,  the  model  will  be  eonstructed  so  that  it  is  universally 
applicable  to  all  elliptical  orbits,  whereas  previous  models  have  been  narrowly  foeused 
on  orbital  regimes  of  interest  (whieh  allows  for  equations  whieh  exhibit  singularities  in 
some  cases  outside  of  study).  As  a  fundamental  premise  of  this  researeh,  it  is  believed 
that  a  validated  single  satellite  model  can  be  replicated  into  several  nearly  identical 
versions  of  itself  and  intereonneeted  in  such  a  way  that  they  ean  all  be  run 
simultaneously.  By  doing  so,  it  is  hoped  that  multi-agent  optimal  control  can  be  shown 
as  nothing  more  than  a  slight  extension  of  well  understood  theories,  and  that,  in  theory,  a 
large  number  of  agents  ean  be  controlled  through  this  process. 
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C.  DIDO  AND  THE  MULTI-SATELLITE  OPTIMIZATION  MODEL 

Although  it  will  be  discussed  very  lihle  through  the  rest  of  this  thesis,  the  DIDO 
optimization  tool  developed  at  the  Naval  Postgraduate  School  by  Dr.  I.  Michael  Ross  and 
Dr.  Fariba  Fahroo  represents  the  enabling  technology  which  has  made  this  research 
possible.  It  is  the  first  and  only  object-oriented  computer  program  for  solving  dynamic 
optimization  programs  [Ross  and  Fahroo  (2002)],  and  it  uses  a  Legendre  pseudospectral 
method  to  perform  this  task.  However,  one  of  the  great  benefits  of  this  tool  is  that  the 
computation  method  is  almost  completely  transparent  to  the  casual  user,  provided  inputs 
are  framed  in  a  structured  manner  understood  by  the  DIDO  interface,  as  explained  in  the 
DIDO  User’s  Manual.  Therefore  its  inner  workings  are  left  to  be  considered  a  black  box 
for  this  research,  and  the  curious  reader  is  invited  to  browse  several  papers  by  Ross  and 
Fahroo  to  discover  more  about  this  tool. 

All  programming  for  this  research  was  performed  in  MATLAB™  release  13 
(version  6.5)  and  release  14  (version  7.0).  All  developed  MATLAB™  codes  supporting 
the  multi-satellite  optimization  model  (hereafter  referred  to  as  “the  model”)  and  a 
common  version  of  the  DIDO  engine  (DIDO  2003)  are  portable  between  the  two  versions 
with  little  apparent  difference  for  this  application,  other  than  a  significant  increase  in 
processing  speed  with  the  later  version.  The  vast  majority  of  this  thesis  will  focus  on 
more  general  mathematics  upon  which  the  model  is  built,  and  will  not  generally  discuss 
actual  code  language.  However,  it  is  assumed  to  be  understood  that  the  mathematical 
generalities  and  specifics  discussed  throughout  this  work  is  fully  implemented  into  an 
operational  MATLAB^^  based  code  structure  wrapped  around  the  DIDO  optimization 
tool. 
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II.  PRINCIPLES  OF  DYNAMIC  OPTIMIZATION 


A.  SOLVING  OPTIMAL  CONTROL  PROBLEMS 

The  driving  principle  used  to  solve  optimal  control  problems  was  first  formalized 
by  the  Soviet  mathematician  Pontryagin  in  the  1950s  [Kopp  (1962)],  and  is  generally 
known  as  the  Minimum  Principlei.  His  work  showed  conclusively  that  optimal  solutions 
to  all  differential  equations  exhibit  several  observable  characteristics,  even  though  the 
solutions  themselves  do  not  readily  announce  themselves  as  optimal.  By  verifying  these 
observables  numerically,  we  can  use  engineering  judgment  to  conclude  that  a  given 
solution  is  in  fact  optimal.  Stated  another  way,  optimal  solutions  must  meet  several 
necessary  conditions  as  proof  of  optimality,  and  by  some  reverse  logic,  if  we  can  observe 
the  necessary  conditions  of  optimality,  we  can  with  fair  confidence  conclude  that  a 
solution  is  optimal.  In  order  to  determine  the  optimal  controls,  several  steps  must  first  be 
taken,  as  detailed  below.2 

1,  Define  the  Performance  Index,  States,  and  Controls 

Bryson  (1999)  defines  optimal  control  as  “the  process  of  determining  control  and 
state  histories  for  a  dynamic  system  over  a  finite  time  period  to  minimize  a  performance 
index.”  In  order  to  solve  for  optimality,  we  must  first  decide  what  index  we  want  to 
minimize  (or  maximize).  The  majority  of  this  research  has  focused  on  minimizing  the 
amount  of  fuel  that  is  consumed  by  the  spacecraft  in  performing  their  combined 
maneuvers.  This  can  also  be  described  as  maximizing  the  end  mass  of  the  satellite.  In 
truth,  maximum  and  minimum  conditions  are  largely  interchangeable  and  are  realizable 
simply  with  a  change  of  sign  (i.e.,  a  minimized  condition  can  also  be  considered  a 
negative  maximized  condition). 

Regardless  of  the  performance  index  chosen,  the  fundamental  problem  is  based  on 
minimizing  (or  maximizing)  the  index  (or  cost).  We  can  generically  state  this  cost  as: 


1  In  truth,  Pontryagin’ s  work  was  known  as  the  Maximum  Principle,  but  for  reasons  to  be  discussed 
the  two  titles  are  largely  interchangeable 

2  The  following  development  largely  follows  the  course  notes  of  Dr.  I.M.  Ross  for  AA3830  - 
Spacecraft  Guidance  and  Control 
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]  =  E{x{t ^),t +  ^  F{x{t),u{t),t^dt 


(2.1) 


Here  E  is  the  end  point  (Mayer)  eost  and  F  is  the  running  (Lagrange)  eost.  As  we 
will  soon  see,  the  researeh  in  this  thesis  foeuses  exelusively  on  Mayer  eost  indiees. 

The  variables  x  and  u  represent  state  and  eontrol  veetors.  This  is  the  standard 
notational  eonvention  that  will  be  used  throughout  this  doeument.  The  state  variables  are 
required  to  eompletely  deseribe  the  eondition  of  the  system  at  any  point  in  time,  and  need 
to  be  ehosen  wisely  and  in  eoneert  with  seleetion  of  the  dynamie  equations.  Care  and 
eonsideration  should  be  given  in  the  ehoiees  of  state  and  eontrol  to  most  aeeurately  model 
the  system  under  study  as  well  as  to  gain  whatever  residual  benefits  a  partieular  set  of 
variables  has  to  offer. 

2,  Develop  the  Dynamic  Equations 

The  next  step  in  modeling  is  to  develop  a  set  of  dynamie  equations  of  motion 
suitable  for  the  task  at  hand.  These  equations  are  time  derivatives  of  the  state  variables 

dx 

denoted  simply  by  x  =  — . 

dt 

A  given  set  of  equations  is  not  a  silver  bullet;  any  set  of  equations  ean  be  used  if 
they  properly  mathematieally  deseribe  the  physies  of  the  problem  at  hand.  Our 
applieation  is  orbital,  and  one  eould  envision  using  equations  based  on  elassieal  elements, 
equinoetial  elements,  Delaunay  elements,  Poineare  elements  [Vallado  (2001)],  or  any 
number  of  representations  one  eould  imagine.  The  ehoiee  of  equations  matters  only  in 
that  they  must  be  eompletely  eapable  of  deseribing  the  dynamies  of  the  seenarios  under 
study. 

The  study  at  hand  is  then  fundamentally  the  determination  of  state  and  eontrol 
pairs  whieh  minimize  the  eost  index,  subjeet  to  the  dynamie  equations: 

Kt)  =  u(t),  t]  (2.2) 

This  determination  neeessarily  involves  two  boundary  values:  an  initial  eondition, 
and  an  end  manifold.  For  this  reason,  this  system  is  termed  a  two  point  boundary  value 
problem. 
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3.  Develop  the  Boundary  Conditions 

A  set  of  starting  values  for  the  state  variables,  x(to)  =  Xo,  must  be  ereated  to 
provide  initial  eonditions  for  the  dynamie  differential  equations.  A  set  of  end  eonditions, 
x{t^)  =  x^,  must  also  be  ereated  to  provide  an  end  manifold  and  to  eomplete  the  two 

point  boundary  value  problem.  Simply  stated,  these  eonditions  determine  the  state  that 
the  model  starts  from,  and  the  state  in  whieh  we  desire  the  model  to  end. 

4.  Develop  the  Path  Constraints 

A  set  of  path  eonstraints,  g,  must  be  developed  to  eonstrain  the  problem  solution 
to  some  less-than-infinite  set  for  eomputational  purposes.  This  set  of  path  eonstraints  ean 
be  a  family  of  possible  variables  whieh  ineludes  eonstraints  on  the  eontrols,  whieh  will  be 
our  primary  eoneem.  Upper  and  lower  path  eonstraints  must  be  set,  where 
g,  <  g{x,u,t)  <  .  These  path  eonstraints  determine  the  maximum  and  minimum  values 

that  the  variables  ean  aehieve,  and  effeetively  bound  the  solution  spaee  of  the  problem. 

5.  Develop  the  Hamiltonian 

The  first  real  step  in  eheeking  optimality  involves  development  of  the 
Hamiltonian  funetion,  H: 

H{x,u,t,X)  =  F{x,u,t)  +  /i{tY  f  {x,u,t)  (2.3) 

Here/l(t)  generieally  represents  Lagrangian  multipliers  assoeiated  with  the 
dynamie  eonstraint  equations  (our  equations  of  motion),  also  known  as  eostates. 

6.  Develop  the  Lagrangian  of  the  Hamiltonian 

In  order  to  find  the  optimal  eontrol  for  this  problem,  at  eaeh  instant  in  time  we 
must  minimize  the  Hamiltonian  with  respeet  to  the  eontrol  veetor  («).  In  order  to  do  this, 
we  next  need  to  form  the  Lagrangian  of  the  Hamiltonian: 

H(x,u,t,A,^)  =  H{x,u,t,A)  +  (/>{tYg{x,u,t)  (2.4) 

Here  generieally  represents  Lagrange  multipliers  assoeiated  with  the  path 

eonstraints  (ineluding  eontrol  duals,  as  we  will  later  see) 
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7,  Apply  Karush-Kuhn-Tucker  (KKT)  Theorem 

By  applying  KKT  theorem  to  the  above  developed  forms,  we  ean  observe  the 
following  eonditions  neeessary  for  optimality: 


du 


<0 
>0 
=  0 

unrestricted 


if 


gi  =g(x,M,0 
g(x,u,t)  =  g^ 
gi  <  g{x,u,t)  <  g^ 

gl=gu 


(2.5) 


(2.6) 


Here  gi  represents  the  lower  bound  on  the  path  eonstraint  (e.g.,  a  minimum 
eontrol  value),  and  gu  represents  the  upper  bound. 


The  Minimum  Prineiple  states  that  an  optimal  solution  must  meet  these  three 
eonditions,  whieh  are  neeessary  but  not  suffieient  in  proving  optimality.  In  observing  the 
above  neeessary  eonditions  then,  we  ean  use  reverse  logie  and  engineering  judgment  to 
eonelude  with  reasonable  eertainty  that  a  partieular  solution  meeting  these  eonditions  is 
optimal.  This  is  the  primary  method  of  proof  of  optimality  used  in  this  researeh,  and  will 
be  a  reeurring  theme  throughout  the  remainder  of  this  thesis. 


B,  OBSERVING  NECESSARY  CONDITIONS  OF  OPTIMALITY 

As  mentioned  above,  an  optimal  solution  must  meet  several  neeessary  eonditions 
aeeording  to  the  Minimum  Prineiple.  The  following  diseussion  will  foeus  on  the  methods 
with  whieh  numerieal  observations  were  performed  using  the  previously  deseribed  tools. 

1.  Feasibility 

Although  this  eondition  does  not  prove  optimality  per  se,  it  is  first  neeessary  to 
show  that  a  given  solution  is  feasible  as  a  potential  optimal  solution.  If  the  solution  does 
not  meet  the  feasibility  test,  there  is  little  point  in  earrying  through  the  rest  of  the 
neeessary  ealeulations  to  determine  optimality.  Beeause  the  equations  of  motion  used  in 
this  thesis  are  ordinary  differential  equations,  the  easiest  way  of  showing  feasibility  is  to 
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independently  propagate  a  solution  from  initial  conditions  to  see  how  closely  the  solution 
matches  model  results.  Since  MATLAB™  provides  the  processing  engine  around  which 
our  model  is  wrapped,  the  on-board  differential  equation  solvers  generally  known  as  the 
ODE  toolset  provide  a  natural  choice  for  this  task.  Although  there  are  a  number  of  these 
tools  for  various  conditions  using  different  methods,  it  was  found  that  ODE_45  works 
sufficiently  well  for  this  application.  By  graphing  the  performance  of  the  model  against 
an  independently  propagated  solution,  it  could  be  very  quickly  determined  visually 
whether  or  not  a  solution  provided  a  feasible  answer.  A  representative  example  of  this 
graphical  technique  is  shown  as  Eigure  1  in  order  to  begin  familiarizing  the  reader  with 
the  process  to  be  used  later  in  the  document. 


Figure  1  Sample  graph  of  solution  behavior  vs.  independent  propagation  (for  six 

controls) 


2,  Behavior  of  the  Hamiltonian 

Eor  optimal  control  problems,  the  Hamiltonian  follows  a  simple  rule  as  dictated 
by  the  Hamiltonian  Evolution  Equation,  where; 


dH  _dH 
dt  dt 


(2.7) 
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For  the  types  of  problem  under  eonsideration  this  value  generally  is  equal  to  zero, 
indieating  that  the  Hamiltonian  maintains  a  eonstant  value  over  the  optimal  time  span 
dH 

(since  - =  0 ).  This  development  will  be  shown  once  specifics  for  this  family  of 

dt 

problems  are  discussed. 

The  actual  value  of  the  constant  value  Hamiltonian  depends  on  the  type  of 
problem  solved.  The  value  is  dependent  on  a  terminal  transversality  condition  called  the 
Hamiltonian  Value  Condition,  which  can  be  stated  as: 

dF 

<2.S, 

Here,  E  is  the  end-point  Mayer  cost,  and  E  is  developed  as  follows: 

E  =  EF^v^e, 

/=! 

where 

=  number  of  equations  defining  the  endpoint  manifold  (2.9) 

e.=x.{tf)-x. 

V-  =  end  state  Lagrange  multiplier 


As  we  will  see  in  a  later  chapter,  the  value  of  the  Hamiltonian  is  generally  zero 
for  minimum  fuel  problems.  Using  this  fact  coupled  with  the  Hamiltonian  Evolution 
equation,  it  can  be  concluded  that  the  Hamiltonian  will  remain  zero  over  the  optimal 
control  period.  For  minimum  time  problems,  this  constant  value  is  -1.  As  with  the 
feasibility  analysis,  a  quick  graphical  method  will  be  used  to  verify  compliance  with  this 
condition,  and  a  representative  graph  sample  is  shown  as  Figure  2. 
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Figure  2  Representative  Hamiltonian  behavior  for  optimal  solution 

3,  Minimized  Lagrangian  of  the  Hamiltonian  with  Respect  to  Controls 

As  discussed  earlier,  it  is  neeessary  that  KKT  equation  (2.5)  must  be  equal  to 
zero  for  optimality.  This  equation  represents  the  stationary  Lagrangian  of  the 
Hamiltonian  with  respeet  to  eontrol  (for  sake  of  simplieity,  this  will  be  alternately 
referred  to  for  the  rest  of  this  thesis  as  the  stationary  Lagrangian).  This  equation  must 
hold  true  for  eaeh  eontrol,  and  therefore  the  number  of  equations  implied  is  equal  to  the 
number  of  eontrols  used.  By  representing  this  eondition  graphieally  for  each  of  the 
eontrols  under  eonsideration,  we  ean  quiekly  get  a  feeling  for  whether  or  not  the 
eondition  meets  our  optimality  test,  namely  that  a  eonstant  zero  should  be  maintained 
throughout  the  time  span.  A  representative  example  of  this  graphical  technique  is  shown 
in  Figure  3. 
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Figure  3  Representative  behavior  of  the  stationary  Lagrangian  of  the 
Hamiltonian  with  respect  to  controls  for  optimal  solution  (six  controls) 

4.  Complementarity  Condition 

As  described  by  equation  (2.6),  a  relationship  exists  between  the  Lagrange 
multipliers  and  the  controls  under  observation.  By  graphically  representing  these 
together,  a  definite  pattern  of  switching  should  occur.  For  example,  when  the  control 
multiplier  is  less  than  zero,  its  associated  control  should  be  at  its  minimum  value. 
Similarly,  when  the  multiplier  is  greater  than  zero,  the  control  should  be  at  its  maximum. 
When  the  multiplier  crosses  through  zero,  the  control  should  demonstrate  a  switch  from 
one  condition  to  the  other,  consistent  with  the  direction  of  the  crossing.  Again,  a 
graphical  technique  will  be  used  for  compliance,  and  a  representative  example  is  shown 
as  Figure  4. 


12 


Some  further  explanation  of  eontrol  switching  behavior  and  model  behavior  is 
warranted  here.  One  of  the  input  parameters  used  with  the  model  is  the  number  of 
discrete  time  nodes  to  be  used  for  a  given  scenario  run.  This  can  be  roughly  equated  as  a 
measure  of  resolution  with  which  the  discrete  solution  can  approximate  a  continuous  time 
result.  For  the  sake  of  processing  speed,  the  number  of  nodes  can  be  set  low.  When 
observing  the  switching  behavior,  this  can  sometimes  have  the  adverse  effect  of  making 
the  control  appear  as  if  it  is  not  achieving  maximum  or  minimum  (bang-bang)  control, 
but  that  the  control  is  throttling  in  some  manner.  When  the  number  of  nodes  for  a  given 
scenario  is  increased,  bang-bang  control  is  more  evidently  displayed. 

Generally,  the  nodal  number  for  the  scenarios  studied  is  set  at  what  is  to  be 
considered  mid-range.  In  some  cases,  bang-bang  control  is  not  immediately  evident, 
however,  when  the  number  of  nodes  is  increased,  bang-bang  control  can  generally  be 
observed. 


5,  Others 

Although  the  above  mentioned  conditions  were  primarily  used  as  verification  of 
optimality  for  modeling,  there  are  other  conditions  which  were  deemed  overly  complex 
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to  prove  for  little  residual  return  for  the  already  complicated  equations  of  motion  under 
consideration  in  this  research.  These  conditions  include  development  of  adjoint 
equations  and  determination  of  further  terminal  transversality  conditions,  which  can  be 
reviewed  in  greater  depth  in  the  text  by  Bryson  and  Ho  (1975). 


C.  SCALING 

One  of  the  peculiarities  of  working  with  optimal  control  codes  is  that  they 
perform  best  when  the  inputs  are  well  scaled  relative  to  one  another  [Ross  and  Fahroo, 
(2002)].  The  inputs  used  should  be  numerically  close  to  each  other  in  an  attempt  to  avoid 
calculations  involving  vastly  different  number  ranges  which  might  throw  answers  outside 
the  bounds  of  computational  precision.  By  scaling  inputs  prior  to  computation,  we  can 
avoid  these  problems  entirely  as  long  as  we  maintain  the  understanding  that  outputs  of 
the  system  will  also  be  in  scaled  variables.  This  can  be  easily  rectilied  by  reversing  the 
scaling  process  at  the  end  of  the  computation  to  achieve  results  in  the  same  units  as  the 
original  inputs. 

The  choice  of  a  particular  set  of  state  variables  can  go  a  long  way  in  assisting  with 
this  scaling  problem.  If  the  elements  of  the  state  vector  are  naturally  well  scaled,  matters 
are  simplified  somewhat,  as  will  be  demonstrated  in  the  following  chapter. 


D,  CONCLUSIONS 

The  basic  theory  of  the  Minimum  Principle  has  been  summarized,  and  the 
following  chapters  will  show  how  this  theory  will  be  applied  to  the  specific  topic  of  this 
research.  The  tenets  of  the  Minimum  Principle  will  be  used  throughout  this  thesis,  and 
form  the  basis  for  all  verifications  of  optimality  that  will  be  discussed. 
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III.  THE  EQUINOCTIAL  ELEMENT  SET 


In  order  to  build  a  model  whieh  is  universally  applieable  to  any  orbit  transfer 
problem,  we  need  to  define  equations  of  motion  whieh  ean  be  applied  to  all  potential 
eonditions.  These  equations  ean  be  defined  in  any  way  that  ultimately  deseribes  orbital 
motion  of  a  body  around  a  eentral  gravitational  foree;  however,  we  want  to  simplify  this 
as  mueh  as  possible  so  that  no  “speeial  oases”  exist  that  fall  outside  the  ohosen  equations 
and  oomplioate  the  algorithms. 


A,  DRAWBACKS  OF  THE  ORBITAL  ELEMENT  SET 

One  of  the  most  oommonly  used  element  sets  in  orbital  meohanios  is  the  olassioal 
orbital  element  (COE)  set.  Although  this  set  is  fairly  easily  understood,  it  has  well 
known  oases  for  whieh  it  exhibits  singularities  [Vallado  (2001)].  In  partioular,  oertain 
elements  beoome  undefined  in  either  perfeotly  oiroular  orbits  or  equatorial  orbits.  For 
zero  eooentrioity  orbits,  the  argument  of  periapsis  has  no  meaning.  For  zero  inclination 
orbits,  the  right  ascension  does  not  exist.  Substitutions  for  these  “special  cases”  must  be 
accounted  for  in  order  to  have  a  universally  applicable  tool.  Although  scaled  to  canonical 
values,  Delaunay  elements  exhibit  similar  difficulties  with  these  special  cases.  While  the 
difficulties  are  not  large  enough  to  discount  these  element  sets  from  use  (workarounds 
can  be  coded  into  the  model),  these  reasons  as  well  as  a  new  challenge  drove  the  early 
decision  to  research  the  use  of  an  alternate  element  set  which  did  not  pose  these 
problems,  namely  the  equinoctial  orbital  element  (EOF)  set. 

B,  DEFINITION  OF  THE  EQUINOCTIAL  ELEMENTS 

An  overview  of  the  equinoctial  elements,  widely  credited  to  Broucke  and  Cefola 
(1972)  for  initial  exploration,  is  provided  in  this  section.  Although  several  conventions 
have  been  used  in  various  published  discussions  of  these  elements,  the  convention  used 
by  Battin  (1999)  has  been  adopted,  and  a  more  detailed  derivation  of  these  elements 
exists  there.  Although  these  elements  can  be  applied  to  parabolic  and  hyperbolic  orbits 
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with  some  variation  [Coverstone  and  Prussing  (2003)],  this  version  of  the  model  has  been 
limited  to  elliptical  orbits  (including  circular  orbits)  only. 

To  begin  describing  the  equinoctial  element  set,  the  reference  frame  will  first  be 
defined.  The  equinoctial  elements  are  best  defined  within  the^w  equinoctial  frame  as 
illustrated  in  Figure  5.  This  frame  can  be  described  as  a  three  rotation  sequence  from  the 
xyz  frame:  a  positive  Q  rotation  about  z  ,  a  positive  i  rotation  about  / ,  and  a  negative  £2 
rotation  about  w . 


Figure  5  Equinoctial  reference  frame 

There  are  six  elements  in  the  equinoctial  element  set,  listed  below  and  depicted  in 
Figure  6: 

•  a  =  semimajor  axis 

•  Pi  =  g-component  of  the  eccentricity  vector 

•  P2  =  f-component  of  the  eccentricity  vector 

•  Qi  =  g-component  of  the  ascending  node  vector 

•  Q2  =  f-component  of  the  ascending  node  vector 
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L  =  true  longitude 


In  order  to  help  understanding,  we  ean  define  these  elements  as  they  relate  to  the 
elassieal  elements:  semimajor  axis,  eoeentrieity,  inelination,  right  aseension  of  aseending 
node,  argument  of  perigee,  and  true  anomaly  {a,e,i,Q.,a),v),  as  diseussed  by  Danielson 
et  al(1995). 

The  semimajor  axis  for  the  EOE  set  is  the  same  as  that  of  the  COE  set,  in  line 
with  standard  eonie  geometry,  and  will  be  assumed  to  require  no  further  explanation. 

The  Pi  and  P2  elements  together  have  a  magnitude  equal  to  the  eeeentrieity  of  the 
orbit  and  form  a  veetor  that  points  to  periapsis  from  the  gravitational  eenter,  where: 


=  esin(n;  +  /Q) 

(3.1) 

P2  =  ecos{co  +  ID.) 

(3.2) 

The  variable  /  as  used  here  is  a  retrograde  faetor,  whieh  is  +1  for  direet  orbits  and 
-1  for  retrograde  orbits.  However,  this  is  ignored  by  many  texts  sinee  the  -1  faetor  is  only 
required  for  true  retrograde  orbits  (i.e.  -  exhibiting  180°  inelination),  whieh  are  very 
rarely  used. 
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The  Qi  and  Q2  elements  together  have  a  magnitude  dependent  on  the  inelination 
and  form  a  veetor  that  points  to  the  aseending  node  from  the  gravitational  eenter,  where: 


Qi  =  tan 


sinQ 


(3.3) 


=  tan 


cosQ 


(3.4) 


The  L  element  defines  position  of  the  satellite  in  the  orbit,  and  is  considered  the 
only  rapidly  changing  variable  in  this  element  set.  It  is  related  to  the  classical  true 
anomaly,  right  ascension,  and  argument  of  perigee  by  the  following: 

L  =  <d  +  v  +  IQ.  (3.5) 


C.  FEATURES  OF  THE  EQUINOCTIAL  ELEMENT  SET 

Equinoctial  elements  are  a  very  good  candidate  for  this  thesis  for  reasons 
discussed  below.  Unfortunately,  they  can  be  mathematically  challenging  to  use  in 
practice,  and  this  perhaps  remains  their  principal  drawback.  However,  through  careful 
and  thorough  coding,  this  problem  can  be  overcome. 

1.  Elimination  of  Singularities 

The  primary  reason  EOEs  were  chosen  for  the  model  is  that  they  do  not  suffer  the 
singularity  problems  exhibited  by  the  COE  set  and  others  in  dealing  with  circular  and 
equatorial  orbits.  Since  (ideally)  equatorial  orbits  and  (ideally)  circular  orbits  play  such  a 
prominent  role  in  the  satellite  systems  that  we  use,  it  was  desired  to  model  using  a  system 
which  could  easily  handle  these  cases.  In  fact,  the  only  orbit  for  which  the  principal 

equinoctial  elements  exhibit  a  singularity  is  direct  retrograde  (i.e.,  -  exhibiting  180° 
inclination),  and  this  can  be  rectified  by  using  a  retrograde  factor  as  discussed 
previous  ly3. 


3  Similarly,  using  the  retrograde  factor  at  -1,  the  only  orbit  where  the  retrograde  elements  exhibit 
singularity  is  for  inclination  of  0°. 
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2,  Scaling 

As  discussed  previously,  sealing  is  a  major  eonsideration  when  performing 
optimal  eontrol  eomputations,  and  using  naturally  well  sealed  state  variables  has  instant 
benefit.  Equinoetial  orbital  elements  are  sueh  variables,  and  with  the  possible  exeeption 
of  the  semimajor  axis,  the  rest  of  the  variables  are  already  well  sealed.  Several  eonstants 
that  are  used  globally  within  the  modeling  eode  also  need  to  be  sealed  in  order  to  ensure 
effieient  eomputation.  The  following  is  a  diseussion  of  how  sealing  is  performed  for  the 
model  used  in  this  researeh. 

First,  four  primary  eonstants  are  defined  whieh  seed  the  sealing  proeess. 
Although  these  values  remain  eonstant  for  a  given  run  of  the  model,  these  values  are 
modeled  as  variables  for  ease  of  modifieation  for  different  runs.  The  eonstants  used  are 
thruster  exit  veloeity  (Ve=Isp*go;  where  speeifie  impulse  (Ep)  and  Earth’s  gravitational 
aeeeleration  (go=9.81*10'  km/see)  are  inputs  not  used  beyond  this  ealeulation),  max 
thrust  (Tmax),  the  gravitational  eonstant  for  an  Earth  eentered  system  (po=398600.4415 
km  /s  ),  and  the  initial  mass  of  the  satellite  (Mq).  These  eonstants  serve  as  the  basis  for 
further  sealing,  and  are  used  in  sealed  form  globally  within  the  model. 

Next,  generie  sealing  faetors  are  developed  to  eonvert  physieal  units  to  sealed 
units.  The  five  prineipal  sealing  faetors  are: 

•  Du  (Distanee  Unit)  =  Radius  of  the  Earth,  6378.1363  km 

•  Mu  (Mass  Unit)  =  Mq 

•  Tu  (Time  Unit)  = 

•  Fu  (Foree  Unit)  = 

Tu 

F)jj 

•  Vu  (Veloeity  Unit)  = - 

Tu 

Using  these  faetors,  sealed  versions  of  the  prineipal  eonstants  are  developed  for 
further  use  throughout  the  model  eode: 
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V 

Scaled  Ve  =  — ^ 

(3.6) 

Vu 

Scaled 

Fu 

(3.7) 

Tu^ 

Scaled  fio  =  //q  — % 

Du 

(3.8) 

M 

Scaled  Mo-  “ 

(3.9) 

Mu 

As  mentioned  above,  the  equinoetial  orbit  elements  are  naturally  suited  to  sealing. 
With  the  exeeption  of  the  semimajor  axis,  whieh  when  provided  in  kilometers  must  be 
divided  by  the  Du  seale  faetor  to  normalize,  the  remainder  of  the  elements  are  well 
sealed.  The  Pi  and  P2  variables  are  eonverted  from  COEs  using  sin  and  eos  funetions  and 
therefore  fall  between  0  and  1  automatically.  The  Qi  and  Q2  variables  are  converted 
using  a  tan  function,  and  although  they  are  well  behaved  for  most  values,  they  exhibit 
difficulty  as  inclination  approaches  180°  (when  1=1)  and  must  be  limited  by  enforcing 
bounds,  as  will  be  discussed  in  the  following  chapter.  The  true  longitude  is  measured  in 
radians  in  orbit  and  therefore  increases  slowly  enough  (2n  per  orbit)  not  to  be  an  issue 
over  relatively  short  periods. 

3,  Mathematical  Complexity 

Unfortunately,  the  benefits  of  EOEs  do  not  come  without  a  price.  Working  with 
the  elements  can  quickly  become  very  complicated  mathematically,  as  evidenced  in  the 
works  of  Betts  (2001),  Kechichian  (1990,  1991),  and  others.  Although  these  sources 
provide  excellent  treatments  of  the  complicated  matrix  mathematics  required  for  classical 
application,  the  use  of  the  tools  involved  in  this  research  have  made  replication  of  most  of 
this  work  unnecessary  outside  of  a  formulation  of  the  equations  of  motion.  This  in  itself 
has  proven  to  be  one  of  the  major  benefits  of  this  research,  and  will  be  covered  in  more 
detail  in  the  following  chapter. 
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D,  TRANSFORMATION  OF  EQUINOCTIAL  ELEMENTS  TO  POSITION 
AND  VELOCITY 

Transformations  of  equinoctial  elements  to  classical  elements  can  be 
accomplished  simply  by  reversing  the  equations  provided  earlier  in  this  chapter. 
However,  for  some  applications,  transformation  directly  to  position  and  velocity  vectors 
is  preferred.  The  following  equations  [Battin,  (1999)]  provide  the  means  of  this 
transformation. 


r  =  r 


COST 

sinT 


0 


V 


-/]  -  sin  L 
P2  +  cos  L 


(3.10) 


(3.11) 


where; 

h  =  angular  momentum  =  “^2^) 

a(\-P^  -P^) 

r  =  radius  from  gravitational  center  = - ^ - - - 

l  +  Py  sin(T)  +  P2  cos(T) 

p  =  semiparameter  =  a(l  -P^  -  P2) 

//g  =  gravitational  constant 


A  rotation  matrix  must  be  used  to  transform  these  resultant  vectors  from  fgw 
coordinates  to  xyz  coordinates; 


R  = 


1 

i+a'+a' 


1-01^ +e/ 
-201 


2002  20 
1  +  01^-02'  -202 
202  l-0^-02' 


(3.12) 


Using  these  together; 


r 

{xyz  frame} 

y 

{xyz  frame} 


=  R*r 

*  {fgw  frame} 

=  R*  V 

'  {fgw  frame} 


(3.13) 
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E,  CONCLUSIONS 

The  judicious  selection  of  state  variables  can  provide  much  benefit  to  a  particular 
problem.  For  this  research,  equinoctical  orbit  elements  have  been  chosen  due  to 
favorable  scaling  properties  and  lack  of  singularities  for  all  orbits.  The  following  chapter 
will  describe  how  these  states  will  be  employed  into  dynamic  equations  completely 
describing  orbital  motion,  and  form  the  foundation  for  all  of  the  research  that  follows. 
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IV.  FORMULATION  OF  THE  PROBLEM 


In  order  to  move  away  from  the  broad  theory  that  has  been  diseussed  thus  far  and 
work  into  the  speeifie  applieations  for  the  problem  under  study,  a  speeifie  optimal  eontrol 
problem  used  throughout  this  researeh  will  be  formulated.  The  following  development 
will  be  performed  for  a  single  agent;  however  as  will  be  later  shown,  the  desired  multi¬ 
satellite  problem  is  a  relatively  simple  extension  of  this  basie  problem. 


A,  STATES  AND  CONTROLS 

As  diseussed  earlier,  the  states  are  the  equinoetial  orbital  elements.  The  state 
veetor  is  symbolieally  expressed  as  follows: 


f  a 


Pr 

P2 


x  =  < 


a 


a 


L 

M 


(4.1) 


Mass  is  required  to  represent  ehanges  as  fuel  is  eonsumed  by  the  satellite’s 
thrusters. 

The  eontrol  variables  for  this  problem  were  deeided  to  be  six  body-fixed  thrusters: 
one  for  eaeh  direetion  in  the  +/-  XYZ  seheme.  This  deeision  was  in  part  based  on 
eonsiderations  of  using  the  L'  eontrol  norm  for  mass  flow  ealeulation,  and  follows 
diseussion  by  Ross  (2004)  regarding  minimum  fuel  eontrollers.  Thus,  the  eontrols  used 
for  this  problem  are: 
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(4.2) 


Here  T  represents  thrust,  r,  t,  and  n  represent  radial,  transverse,  and  normal  direetions, 
and  the  1  and  2  subseripts  represent  positive  and  negative  directions  respectively.  Thus, 
r^i,7^i,r„i  >  0  while  <  0  .  A  diagram  of  this  convention  is  shown  in  Figure  7. 


Tm 

Figure  7  Satellite  thrust  convention  diagram.  Note,  positive  normal  thrust  T„i  not 

shown 


B,  COST 

The  main  goal  of  this  research  is  to  minimize  the  amount  of  fuel  consumed  during 
commanded  maneuvers.  As  discussed  previously,  this  can  also  be  stated  as  maximizing 
the  final  mass  of  the  vehicle.  If  we  consider  the  Bolza  cost  functional,  then  only  the  end¬ 
point  Mayer  cost  E  needs  to  be  considered,  so  that: 

J  =  =  (4.3) 

Here  the  state  variable  M  is  the  mass  of  the  satellite  at  any  point  in  time,  and  Mf  is  the 
final  mass  of  the  vehicle,  the  variable  to  be  maximized. 
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C.  DYNAMIC  EQUATIONS  OF  MOTION 

The  maximization  of  the  performanee  index  A//^must  be  accomplished  subject  to 
dynamic  equations  governing  orbital  motion.  These  equations  are  based  on  the 
derivatives  of  the  state  variables,  and  can  be  represented  simply  as; 


Pi 


x  =  < 


a 


\  =  f{x,u) 


4 


L 

M 


(4.4) 


The  right  hand  sides  of  seven  the  differential  equations,  f(x,u),  are  intricately 
coupled.  As  stated  previously,  the  derivation  of  the  full  equations  is  based  on  Battin’s 
text  (1999)  coupled  with  the  chosen  controls.  The  full  equations  of  motion  can  be  written 
as  follows: 


d  =  2 


V  "  7 


(P,sin(i)-i^cos(i))| 


P 
+  1  - 


T  +T 

-'<1  ^-'<2 

M 


(4.5) 


^=1- 


.4 

r  J 


T  +T 

cos(T)|  — - — 

'  M 


+ 


1  +  1^ 


P^{Q^coi{L)-Q^?,m{L)) 


T  +T 
M 


r  JJ 

)) 


sin(Z) 


\  /'  rp  p  \\ 

-'<1  ^  -*<2 

M 


- 


r  J 

f 


sin(Z) 


T  +T 

-'rl  ^  -'r2 

M 


P2  + 


1  +  1^ 


i^(giCOs(Z)-g2sin(Z)) 


T  +T 

-'nl  ^  -'ll 

M 


r  JJ 

\\ 

2^ 

)) 


cos(T) 


T  +T 

-'ll  ^  -'<2 

M 


r  r  \ 


01=^— J(1 +  21+02  )sin(T) 


T  +T 
M 


02=|^J(l  +  ei  +02  )cos(T) 


T  +T 

^n\  ^  -^nl 

M 


(4.6) 


(4.7) 


(4.8) 

(4.9) 
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(4.10) 


^h4 

r  J 


+ 


rh 


(asin(Z)-acos(Z))f^44^^ 

\  Tkz 


M 


iT,r-T.^)  +  {T„-T,,)  +  {T„,-T„,) 


where; 

h  =  angular  momentum  =  “^2^) 

_/52  _p2x 

r  =  radius  from  gravitational  center  ~  ^ 


1  +  sm(Z)  +  cos(Z) 


p  =  semiparameter  =  a{\- P^  - P^) 
//q  =  gravitational  constant 
Vg  =  exit  velocity  of  thruster 


(4.11) 


D,  EVENTS 

Since  this  is  a  two  point  boundary  value  problem,  the  only  events  required  are  the 
initial  and  final  boundary  conditions  for  the  state  variables.  Further,  since  we  are  looking 
to  maximize  final  mass,  Mf  can  be  left  undefined.  Thus,  our  relevant  boundary  events 
are; 

a(to)  =  a. 

^2(^0)  ~  1^2i 
2i(^o)  ~  Q\i 
62(^0)  =  62, 

L{t,)  =  L, 

M{Q  =  M, 

Outside  of  some  path  constraints  and  defining  constants,  these  are  the  only  inputs 
required  by  the  model  to  determine  optimal  control. 


a(t^)  =  Uj- 

^  (^/  )  “  ^/ 
P^itf)  =  P2J 

Qiitf)  =  Qif 

Q2i(f)  =  Q2f 
L{t  =  L  ^ 


E,  STATE  AND  CONTROL  BOUNDS 

Some  boundaries  need  to  be  enforced  on  states  and  controls  in  order  to  define  a 
subset  of  possible  solutions  inside  of  which  we  will  find  a  locally  optimum  one.  This 
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helps  to  constrain  a  potentially  mathematically  infinite  range  of  possible  solutions  to 
some  reasonable  (but  not  overconstrained)  subset.  Each  of  the  set  boundaries  bears  some 
discussion,  but  state  bounds  have  been  set  as  follows.  Note  that  all  values  shown  are 
scaled  as  previously  discussed. 

1,  State  Bounds 

The  bounds  on  semiparameter  are  straightforward,  and  are  simply  1  Earth  radius 
as  the  lower  bound  (to  avoid  “Earth  intercept”  orbits),  and  an  arbitrarily  large  number 
(not  infinity  for  computational  reasons)  as  the  upper  bound. 

1<  a  <10000  (4.12) 


The  bounds  on  Pi  and  P2  are  set  in  order  to  constrain  the  possible  set  to  elliptical 
orbits  only,  since  the  two  parameters  are  based  upon  the  orbit’s  eccentricity 

j  •  This  relationship  is  also  defined  as  a  path  constraint  in  order  to 

maintain  ellipticity  within  the  transfer  (i.e.,  0  <  e  <  1 ).  In  reality,  both  Pi  and  P2  can  not 
be  at  the  maximum  (or  minimum)  value  at  the  same  time;  however,  the  bound  is  set  at 
less  than  1  so  that  a  parabolic  orbit  (i.e.,  where  e=l)  is  not  achieved  if  one  value  equals  1 
while  the  other  is  0. 


-.999<P,<.999 

-.999<P^<.999 


(4.13) 


The  bounds  on  Qi  and  Q2  are  again  simply  set  arbitrarily  high.  Since  the 
conversion  to  the  Q  variables  is  reliant  on  a  function  involving  tan(i/2),  most  values  are 
within  reasonable  scaled  ranges.  However,  as  i  approaches  180°,  these  values  begin  to 
grow  towards  infinity,  so  must  be  given  room  to  run. 


-10000  <gi<  10000 
-10000  <02  ^10000 


(4.14) 
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The  bounds  on  true  longitude  are  set  between  zero  and  an  arbitrarily  high  upper 
value  whieh  is  set  as  a  eonstant  within  the  model  eode. 

0<T<  10000  (4.15) 


The  bounds  on  mass  are  set  between  an  arbitrary  non-zero  minimum,  with  the 
initial  mass  of  the  satellite  used  for  the  upper  bound.  The  minimum  ean  be  set  as  a  way 
of  defining  the  usuable  propellant  portion  of  the  satellite  mass  and  must  be  greater  than 
zero  (sinee  equations  (4.5)  -  (4.10)  are  singular  for  M=0). 

(4.16) 


1.  Control  Bounds 

The  bounds  levied  on  the  eontrol  variables  are  fairly  straightforward.  A 
maximum  thrust  is  defined  as  a  eonstant  based  on  the  individual  thruster  being  modeled. 
The  positive  and  negative  direetions  for  each  dimension  are  then  just  set  by  convention  as 
being  between  0  and  the  max  thrust  for  that  direction.  The  upper  and  lower  bounds  are 
determined  by  which  direction  that  thruster  operates,  as  follows; 


o<7;,<7; 

o<7;,<7; 

o<7;i<7; 

-T^ax  ^  7;2  ^  0 

-T  <  T,  <  0 
-T  <  r ,  <  0 

max  «2 


(4.17) 


F.  PATH  CONSTRAINTS 

For  the  single  agent  model,  only  one  path  constraint  has  been  levied  on  the 
system.  As  has  been  previously  mentioned,  a  path  constraint  was  created  in  order  to 
ensure  that  orbits  remained  elliptical  at  all  times,  including  transfer.  This  was  done  to 
keep  parabolic  and  hyperbolic  trajectories  from  being  considered  during  transfer,  and  to 
assist  computational  efficiency.  Since  this  model  was  built  to  only  consider  the 

mathematics  of  elliptical  orbits,  this  constraint  serves  to  keep  solutions  within  an 
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understood  regime  and  avoid  wandering  into  territory  for  which  the  system  has  not  been 
suitably  tested^.  The  path  constraint  simply  stated  is; 

<1  (4.18) 

As  we  will  show  later,  additional  path  constraints  must  be  added  once  we  consider 
multi-agent  models  in  order  to  ensure  that  collision  avoidance  between  agents  can  be 
maintained,  as  well  as  to  provide  relational  operating  windows  between  agents. 


G.  DEVELOPMENT  OF  THE  HAMILTONIAN 

The  Hamiltonian  of  the  subject  system  of  equations  can  be  represented  in  a  fairly 
straightforward  manner  (avoiding  full  expansion  for  the  time  being); 

H  =  -Mf+\(d)  +  A,^(P,)  +  A,^(P,)  +  AQ^(a)  +  \SQ2)  +  \U)  +  MM)  (4.19) 


Here  A,  represents  Lagrangian  multipliers  of  the  state  variables,  as  previously 
discussed.  Since  none  of  the  dynamic  equations  are  explicitly  dependent  on  time,  and 
following  the  previously  discussed  Hamiltonian  Evolution  equation  (2.7),  we  find; 


dH  _  dH 
dt  dt 


(4.20) 


The  natural  conclusion  to  this  statement  is  that  the  Hamiltonian  must  exhibit 
constant  behavior  at  all  points  in  time  during  the  optimal  span,  since// =  0.  Further, 
following  the  discussion  of  the  Hamiltonian  Value  condition  (2.8),  we  find; 


4  However,  by  altering  the  semimajor  axis  state  to  semiparameter  (p=a(l-e^)),  it  is  believed  that  this 
model  ean  be  developed  to  work  with  hyperbolie  orbits  as  well 
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H{tf)  = 


BE 


E=E  +  i/,e,  =  -Mf  +  +  Vp^  ep^ 


"’’'^62^62 


e^  =  a{t^)-a 
=Piitf)-Pi 
Bp^  =  P^itf)  —  P2 
eg,  =Qi{tf)-Qi 

“  Ql  (^/  )  “  02 

=L(t^) 


So: 


H{tf)  =  -—  =  Q 

^  dt. 


(4.21) 


Since  the  Hamiltonian  is  zero  at  the  end  state,  and  we  have  shown  that  the 
Hamiltonian  must  remain  constant  throughout  the  optimal  control  span,  we  must 
conelude  that  a  constant  zero  must  be  maintained  throughout  the  eontrol  span  for  a  given 
solution  to  this  problem  to  be  considered  optimal. 


H.  DEVELOPMENT  OF  THE  LAGRANGIAN  OF  THE  HAMILTONIAN 

Using  the  above  results,  we  ean  now  express  the  Lagrangian  of  the  Hamiltonian 
as  follows: 

H  =  -M  ^  +  X^{d)  +  Xp^{P^)  +  Xp^  (P2 )  +  (2i )  +  Xq^  (Q2  )  +  Xp  (Z)  +  (M) 

Pr\  E-T,^  Ptl  P/il  Prl  Et,2  Ptl'^  /^r„2  Zi2 


Here  fj.  represents  Lagrangian  multipliers  on  the  eontrol  variables,  as  previously 
diseussed.  Aeeording  to  KKT  theory,  the  Lagrangian  of  the  Hamiltonian  is  stationary 
with  respeet  to  the  optimal  eontrol  at  all  points  in  time  over  the  optimized  interval,  as 
previously  shown  in  equation  (2.5)  ,  so  that: 


du 


(2.5) 
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If  we  expand  this  relationship  to  its  full  glory,  we  get  the  following  set  of 
equations  which  represent  the  stationarity  of  the  Lagrangian  of  the  Hamiltonian  with 
respect  to  the  controls,  and  which  will  be  used  for  further  proof  of  optimality: 


=  ^  —  (P2sinI-/^cosI)  — ^r^cosl 

Ml  h  \  Mlh 

+  —  — sinZ  +  =0 

M  h  '' 


(4.23) 


=  sinI-/^cosI)]-^r^cosI 

Ml  h  \  Mlh 

+  — ^  — sinZ  + ju^  =0 

M  h  v„ 


(4.24) 


dT^,  Ml  h  yr)\  MU 


Ur  (  p\  ]  A 

-  P,+  1+^  cosL  — ^  +  ju^  =0 
M  /z  M  r  w.  •' 


(4.25) 
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again  where; 

h  =  angular  momentum  =  -  Pi) 

r  =  radius  from  gravitational  center  =  - ^ - - - 

1  +  sin(Z)  +  cos(Z) 

p  =  semiparameter  =  a{\- - P^) 

Pq  =  gravitational  constant 
Vg  =  exit  velocity  of  thruster 


It  should  be  noted  that  these  equations  look  almost  identical  for  any  given  +/- 

X 

thruster  pair,  the  only  differences  being  the  sign  on  the  term,  as  well  as  the  differing 
values  of  the  individual  Lagrangian  control  multipliers. 


I.  COMPLEMENTARITY  CONDITION 

At  each  instant  of  time,  the  multiplier-control  pair  must  satisfy  the  KKT 
complementarity  condition,  which  defines  the  switching  structure  of  the  optimal  control: 


<0 

>0 


[unrestricted 


U:  =  U„ 


if 


u,  =  u. 


U-i  <U<  U-^ 


(4.29) 


As  an  example  of  this,  and  using  previously  discussed  values  for  one  of  the 
controls  (Tri),  we  can  observe: 


Pt  , 


<0  r,!  =  0 

=  0  0<T,<T^ 

unrestricted  0  = 


(4.30) 


By  observing  behavior  following  this  discussion,  we  can  use  this  as  a  third  data 
point  in  validating  optimal  control. 
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J.  CONCLUSIONS 

Having  now  developed  the  speeific  equations,  eonditions,  and  tools  for  aehieving 
and  proving  optimal  control,  we  can  now  focus  on  validating  that  these  building  blocks 
and  the  model  built  upon  them  works  according  to  the  Minimum  Principle. 
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V.  MODEL  VALIDATION 


Having  established  the  basie  methods  and  proeesses  of  solving  optimal  eontrol 
problems,  this  seetion  will  demonstrate  how  the  formulations  work  as  expeeted. 

Onee  eoded,  the  model  was  validated  by  running  well  doeumented  simple  optimal 
eontrol  seenarios  and  observing  the  results  for  agreement.  The  rest  of  this  ehapter 
diseusses  the  efforts  performed  and  their  results  in  eonfirming  optimal  eontrol  behavior 
by  the  model. 

A,  MODEL  NOTES 

As  has  been  alluded  to  previously,  the  MATLAB^m  based  model  aeeepts  inputs  in 
the  form  of  Keplerian  orbital  elements,  then  automatieahy  performs  the  eonversion  to 
equinoetial  elements,  ineluding  automatie  sealing  for  eomputational  efficieney.  This  was 
deemed  the  most  appropriate  and  easily  understood  method  of  input,  sinee  elassieal 
elements  are  eonsidered  somewhat  of  a  basie  standard  among  astrodynamieists. 
However,  onee  the  eonversion  is  made  to  equinoetial  elements  by  the  model,  ah  internal 
ealeulations  are  in  sealed  equinoetial  variables  for  purposes  of  singularity  avoidanee. 

For  the  purposes  of  ah  of  the  seenarios  run  used  in  this  thesis,  a  uniform  set  of 
constants  was  used,  and  they  are: 

/,,=320  s 

go  =0.00981  km/s" 

V.  =4/go=3.1392km/s 

t_=ion 

//o  =398600.4415  kmV^' 

Mo  =  3000  kg 

During  the  following  scenarios,  no  judgment  was  made  as  to  how  much  fuel  was 
consumed  by  the  spacecraft,  only  proving  that  it  was  the  optimal  amount  required  for  the 
requested  transfer.  This  is  to  say  there  was  no  constraint  placed  on  using  for  example 
90%  of  a  vehicle’s  mass  in  fuel  to  perform  a  given  transfer.  For  the  purposes  of  this 
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research,  this  type  of  constraint  was  not  imposed,  but  would  be  the  job  of  mission 
planners  to  determine  reasonable  fuel-to-payload  ratios. 

Scenario  run  speed  was  not  a  major  focus  of  this  effort,  outside  of  the  desire  to 
generally  speed  up  and  automate  the  process  of  determining  fuel  optimal  maneuvers. 
This  is  not  considered  an  issue,  since  current  execution  times  on  the  order  of  man-days 
are  practiced,  and  this  research  has  produced  run  times  of  minutes  or  less.  Approximate 
times  that  individual  runs  have  taken  will  be  noted  in  the  following  sections.  However, 
these  times  are  to  be  taken  with  some  suspicion  since  they  are  platform  and  software 
version  specific,  and  some  overhead  variance  is  incurred  by  running  on  a  network  based 
Windows  operating  system.  For  purposes  of  general  discussion,  all  scenario  runs  were 
performed  on  a  3Ghz  Intel®  Pentium®  4  notebook  computer  using  MATLAB^m  7.0 
running  DIDO. 


B,  HOHMANN  TRANSFER  SCENARIO 

For  the  first  validation  scenario,  a  classic  Hohmann  transfer  was  performed.  First 
suggested  by  Walter  Hohmann  in  1925,  the  Hohmann  maneuver  is  mathematically 
accepted  as  the  most  fuel  efficient  transfer  between  two  circular  coplanar  orbits, 
involving  two  tangential  impulsive  burns  to  start  and  stop  the  transfer,  and  is  well 
documented  by  Vallado  (1999),  Chobotov  (2002),  Bate,  Mueller,  and  White  (1971)  and 
most  introductory  texts  on  orbital  dynamics.  Although  Hohmann’ s  work  only  considered 
circular  orbits,  further  work  by  Lawden  (1952)  validated  that  tangential  bums  are  also 
optimal  for  some  elliptic  scenarios. 

This  orbital  transfer  was  chosen  as  a  classic  example  of  a  known  two-dimensional 
orbit  transfer  with  a  known  optimal  result.  Performance  compliant  with  this  known 
behavior  will  serve  to  validate  this  model  for  a  two-dimensional  coplanar  transfer,  the 
first  step  in  validating  an  overall  three-dimensional  model. 

It  should  be  noted  that  tme  Hohmann  behavior  is  not  expected  from  the 
simulation  since  the  Hohmann  transfer  involves  impulsive  burns  at  the  tangent  points. 
Since  this  model  involves  discrete  finite  burns,  a  similar  behavior  is  expected,  however 
identical  performance  is  not. 
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1,  Hohmann  Transfer  Input  Parameters 

The  following  Keplerian  inputs  were  provided  to  the  model  for  evaluation: 


a.  =  7015  km 

Qj-  =  9567  km 

o 

o 

II 

o 

o 

II 

■ 

i  =  35  *  rad 

'  180 

if  =35*-^  rad 
^  180 

Q  =30*-^  rad 

180 

Q  =  30*-^  rad 
^  180 

(o  =  45  *  rad 

'  180 

CL>f  =45*  rad 

180 

V.  =  0  *  rad 

'  180 

Vf  =235*-^  rad 
^  180 

The  number  N  of  diserete  nodes  used  by  the  model  for  this  seenario  is  80.  This 
number  is  eonsidered  to  be  reasonably  high,  and  impaets  the  eomputational  run  speed  of 
the  seenario.  However  it  offers  the  benefit  of  inereased  resolution  for  illustrative 
purposes.  Computationally  speaking,  the  model  independently  solves  equations 
involving  the  number  of  state  plus  eontrol  variables  times  the  number  of  nodes  aeross  the 
optimal  time  span: 

#  variables  =  (x  +  m)  *  (5.3) 

For  this  example,  this  equates  to  a  (7+6)*80=1040  variable  problem  being 
optimized  for  each  instant  of  time. 

2,  Hohmann  Transfer  Performance 

Given  only  the  above  end  point  conditions,  and  subject  to  the  cost  and  dynamic 
equations  discussed  earlier,  the  following  performance  was  recorded.  Graphical 
representations  of  the  state  and  control  histories  are  shown  as  Figure  8.  Mass  flow  rate  is 
shown  as  Figure  9.  A  pictorial  plot  of  orbital  performance  is  shown  as  Figure  10.  It 
should  be  noted  that  thrusts  occur  only  in  the  positive  transverse  thruster,  as  expected, 
and  that  the  thrust  points  all  fall  in  the  vicinity  of  the  nominal  tangent  points  for  this 
maneuver.  An  80  node  run  took  approximately  124  seconds  to  execute. 
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Figure  8  Hohmann  transfer  state  and  control  histories 


Figure  9  Hohmann  transfer  mass  flow 
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X 

Figure  10  Hohmann  transfer. 


It  is  obvious  by  observation  that  this  performance  does  not  appear  to  be  a  strict 
two  burn  impulsive  maneuver,  but  as  expected,  a  finite  bum  over  several  discrete  points 
closely  simulating  this  behavior.  Since  limits  have  been  placed  on  the  maximum  amount 
of  thmst  available  within  the  model  representative  of  a  reasonably  sized  chemical 
thmster,  the  thmster  is  never  allowed  to  reach  the  force  required  to  accelerate  the  vehicle 
by  the  amount  needed  for  an  impulsive  maneuver.  These  differences  between  theory  and 
a  more  realistic  model  are  accentuated  by  the  fact  that  the  model  used  is  based  upon 
discrete  time  points  rather  than  a  time  continuum,  which  tends  to  magnify  differences  in 
the  thmsting  behavior,  particularly  at  low  discrete  nodal  resolution.  As  the  number  of 
nodes  is  increased,  the  behavior  becomes  more  like  the  theoretical  solution.  This 
seeming  discrepancy  is  one  of  the  sacrifices  made  for  computational  efficiency  in  the 
Legendre  pseudospectral  method,  but  classical  behavior  can  be  extrapolated  from  these 
approximations,  and  indeed  recognized  as  the  number  of  nodes  increases.  This  will  be 
seen  as  a  continuous  theme  through  the  performance  of  the  scenario  runs  through  this 
research,  but  is  no  cause  for  alarm. 
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3. 


Hohmann  Transfer  Model  Optimal  Behavior 


a.  Hohmann  Transfer  Feasibility  Analysis 

Feasibility  was  shown  by  independently  propagating  the  initial  conditions 
through  the  equations  of  motion  using  interpolated  controls.  A  graphical  representation 
of  the  result  is  shown  as  Figure  11,  and  shows  very  good  behavior.  Note  that  the 
propagated  solution  perfectly  overlays  the  model  solution.  It  can  be  concluded  that  the 
result  is  within  the  feasible  set  of  optimal  solutions. 


b.  Hohmann  Transfer  Hamiltonian  Behavior 

Graphical  representation  of  the  Hamiltonian  during  the  performance 
period  is  shown  as  Figure  12.  In  line  with  our  previous  discussion,  it  can  be  observed 
that  the  Hamiltonian  exhibits  a  near  constant  zero  value  for  the  majority  of  the  run.  By 
this  test  condition,  it  can  be  concluded  that  optimal  control  is  being  performed. 
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Figure  12  Hohmann  transfer  Hamiltonian  behavior 


c.  Hohmann  Transfer  Lagrangian  Behavior 

Graphical  representation  of  the  stationarity  of  the  Lagrangian  of  the 
Hamiltonian  (with  respect  to  controls)  during  the  performance  period  is  shown  as  Figure 
13  for  all  six  controls.  In  line  with  our  previous  discussion,  it  can  be  observed  that  the 
Lagrangian  exhibits  a  near  zero  value  for  the  majority  of  the  run.  By  this  test  condition, 
it  can  be  concluded  that  optimal  control  is  being  performed. 
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Hamiltonian  Minimization 


Figure  13  Hohmann  transfer  minimized  Lagrangian  behavior 

d,  Hohmann  Transfer  KKT  Complementarity  Condition 

Graphical  representation  of  the  switching  structure  for  each  of  the 
eontrol/dual  pairs  is  shown  as  Figure  14,  with  greater  detail  visible  in  Figure  15.  This 
switching  generally  follows  the  criteria  described  by  equation  (4.29),  and  therefore  the 
conclusion  can  be  made  that  optimal  eontrol  is  being  performed. 

Note  in  Figure  15  that  the  switeh  seems  to  be  commanded  before  the 
control  dual  crosses  zero.  It  is  also  observed  that  the  value  of  the  costate  remains  very 
close  to  zero  during  the  thrusting  period.  This  close  zero  approach  by  the  eostate  begins 
to  very  closely  resemble  controFdual  behavior  reported  by  Lawden  (1963)  for  an 
impulsive  maneuver,  particularly  since  a  prolonged  period  of  maximum  thrust  is  not 
observed.  In  these  plots,  the  control  does  not  appear  to  reach  maximum  value  at  the 
second  peak  as  might  be  expected;  however  the  control  dual  is  effectively  zero  at  this 
point  so  KKT  theorem  still  holds.  It  is  believed  that  this  behavior  is  due  to  internal 
tolerance  settings  perhaps  coupled  with  some  scaling  issues,  whieh  is  not  a  cause  for 
immediate  concern  but  which  will  be  addressed  further  in  a  later  section  on  future  work. 
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Figure  14  Switching  structure  for  Hohmann  transfer 


Detailed  Complementarity  Switching  Structure 


Figure  15  Detailed  switching  structure  for  positive  transverse  thrust,  Tti,  including 

expanded  blowup  (bottom).  Behavior  approaches  that  reported  by  Lawden. 
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4,  Hohmann  Transfer  Scenario  Conclusions 

Based  on  the  performanee  of  the  model  and  our  observed  eomplianee  with  several 
of  the  neeessary  eonditions  of  optimality,  we  ean  reasonably  eonelude  that  the  model 
exhibits  optimal  eontrol  for  a  two-dimensional  (i.e.-  eoplanar)  orbit  transfer.  As  a  non- 
numerieal  eheek,  the  performanee  of  the  model  also  meets  our  expeetations  on  how  we 
believe  a  Hohmann  transfer  should  be  performed.  As  an  added  benefit,  this  seenario  has 
verified  that  eireular  orbits  (i.e.  -  zero  eoeentrieity  orbits)  do  not  pose  a  singularity 
problem  to  the  model.  This  is  the  first  step  in  validating  overall  three-dimensional 
eomplianee  with  the  mathematies  and  physies  involved  in  Keplerian  orbits.  It  also 
validates  that  our  derived  equinoetial  based  equations  of  motion  appear  to  be  working  as 
expeeted.  However,  we  must  perform  further  validation  to  ensure  the  system  works  for  a 
three-dimensional  ease. 

C.  INCLINATION  CHANGE  SCENARIO 

Sinee  the  model  has  been  validated  for  a  two  dimensional  seenario,  the  next  step 
is  to  prove  that  it  works  for  a  three  dimensional  seenario.  For  this  reason,  a  simple 
inelination  ehange  was  ehosen  to  be  modeled,  sinee  we  know  that  the  mathematieally 
optimal  control  solution  occurs  at  the  ascending  and/or  descending  nodes.  Additionally, 
an  eccentricity  was  added  for  the  initial  and  final  orbits  to  introduce  a  new  element  into 
the  puzzle.  The  presentation  of  this  scenario  will  follow  much  like  the  presentation  of  the 
previous  scenario,  so  some  of  the  explanatory  details  of  the  methodology  will  be  omitted 
as  they  should  now  be  familiar  to  the  reader. 

1.  Inclination  Change  Parameters 

The  following  COE  inputs  were  provided  to  the  model  for  evaluation: 
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a,  =  9567  km 
e,  =  0.5 

i  =  35  *  rad 

'  180 

Q  =10*-^  rad 
‘  180 

cl>-  =  225  *  rad 
'  180 

K  =  0  *  rad 

'  180 


a f  =  9567  km 
=  0.5 

if  =45*-^  rad 
^  180 

Qf  =10*-^  rad 
^  180 


=  225*  ^ 

rad 

180 

=  359*  ^ 

rad 

180 

(5.4) 


Again,  an  80  node  solution  was  requested. 

2.  Inclination  Change  Performance 

Given  only  the  above  end  point  conditions,  and  subject  to  the  cost  and  dynamic 
equations,  the  following  performance  was  recorded.  Graphical  representations  of  the 
state  and  control  histories  are  shown  as  Figure  16.  Mass  flow  rate  is  shown  as  Figure  17. 
A  pictorial  plot  of  orbital  performance  is  shown  as  Figure  18.  It  should  be  noted  that 
thrusts  occur  only  in  the  normal  direction  in  the  vicinity  of  the  ascending  node,  as 
expected.  An  80  node  run  took  approximately  271  seconds. 

Again,  it  can  be  observed  that  this  scenario  does  not  behave  in  impulsive  burns  at 
the  nodes  as  predicted  by  theory,  but  rather  finite  burn  limitations  still  hold  true  for  the 
model.  The  locations  of  the  bums  as  expected  by  optimal  behavior  is  unmistakable,  and 
the  previous  discussion  of  bang-bang  control  is  still  valid,  indicating  a  period  of 
maximum  thmst  at  the  nodal  point  results  in  optimal  control. 
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Figure  16  State  and  control  histories  for  inclination  change 


Figure  1 7  Mass  flow  for  inclination  change 
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Figure  18  Inclination  change 


3,  Inclination  Change  Optimal  Behavior 
a.  Inclination  Change  Feasibility 

Feasibility  was  again  shown  by  independently  propagating  the  initial 
conditions  through  the  equations  of  motion  using  interpolated  controls.  A  graphical 
representation  of  the  result  is  shown  as  Figure  19  and  shows  very  good  behavior.  We  can 
conclude  that  the  result  is  within  the  feasible  set  of  solutions. 
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b.  Inclination  Change  Hamiltonian  Behavior 

Graphical  representation  of  the  Hamiltonian  during  the  performance 
period  is  shown  as  Figure  20.  In  line  with  our  previous  discussion,  we  observe  that  the 
Hamiltonian  exhibits  a  near  constant  zero  value  for  the  majority  of  the  run.  By  this  test 
condition,  we  can  conclude  that  optimal  control  is  being  performed. 


c.  Inclination  Change  Lagrangian  Behavior 

Graphieal  representation  of  the  stationary  Lagrangian  during  the 
performance  period  is  shown  as  Figure  21  for  all  six  controls.  In  line  with  our  previous 
discussion,  we  observe  that  the  stationary  Lagrangian  exhibits  a  near  zero  value  for  the 
majority  of  the  run.  By  this  test  eondition,  we  can  conclude  that  optimal  control  is  being 
performed. 
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Figure  20  Hamiltonian  behavior  for  inclination  change 


Hamiltonian  Minimization 


Figure  21  Inclination  change  minimized  Lagrangian  behavior 


d.  Inclination  Change  KKT  Complementarity  Condition 

Graphical  representation  of  the  switching  structure  for  each  of  the 
control/dual  pairs  is  shown  as  Figure  22  and  detailed  in  Figure  23.  This  switching 
follows  the  criteria  described  by  equation  (4.29),  and  therefore  we  can  conclude  that 
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optimal  control  is  being  performed.  It  can  be  observed  that  the  only  control  being 
exercised  is  by  the  -n  thruster,  Tn2.  It  can  also  be  observed  in  Figure  23  that  this  behavior 
again  appears  to  very  elosely  approach  Lawden’s  (1963)  dual/control  behavior  for  an 
impulsive  bum.  Again,  manifestations  of  the  sealing/tolerance  behavior  observed  in  the 
Hohmann  scenario  can  be  noted,  although  in  this  seenario  we  do  observe  the  occurrence 
of  maximum  thmst  (0.34  scaled  units). 
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Figure  23  Switching  detail  for  thruster  Tni,  with  blow-up 


4.  Inclination  Change  Scenario  Conclusions 

Based  on  the  performance  of  the  model  and  our  observed  compliance  with  several 
of  the  necessary  conditions  of  optimality,  we  can  reasonably  conclude  that  the  model 
exhibits  optimal  control  for  a  three-dimensional  orbit  transfer.  As  a  non-numerical 
check,  the  performance  of  the  model  also  meets  our  expectations  on  how  we  believe  an 
optimal  inclination  change  should  be  performed.  At  first  blush,  this  seems  to  validate 
that  our  derived  equinoctial  based  equations  of  motion  appear  to  be  working  as  expected. 
However,  we  researched  further  scenarios  to  try  to  get  a  better  feeling  that  this  was  in  fact 
the  case. 

D.  SEMIMAJOR  AXIS  CHANGE  SCENARIO 

Although  two-  and  three-dimensional  model  behavior  had  been  proven,  the  next 
step  was  to  get  further  data  points  on  model  performance  in  known  optimal  problems.  A 
semimajor  axis  change  maneuver  was  attempted,  as  frrst  proposed  by  Lawden  (1962). 

1,  Semimajor  Axis  Change  Parameters 

The  following  COE  inputs  were  provided  to  the  model  for  evaluation; 
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a,  =  9567  km 
e,  =  0.5 

i  =  35  *  rad 

'  180 

Q  =10*-^  rad 
‘  180 

cl>-  =  225  *  rad 
'  180 

K  =  0  *  rad 

'  180 


a j  =  9567  km 
=  0.5 

if  =35*-^  rad 
^  180 

Qn  =  10*-^  rad 
^  180 


=  315*  ^ 

rad 

180 

=  359*  ^ 

rad 

180 

(5.5) 


2.  Semimajor  Axis  Change  Performance 

Given  only  the  above  end  point  eonditions,  and  subjeet  to  the  eost  and  dynamie 
equations,  the  following  performanee  was  reeorded.  Graphieal  representations  of  the 
state  and  eontrol  histories  are  shown  as  Figure  24.  Mass  flow  rate  is  shown  as  Figure  25. 
A  pietorial  plot  of  orbital  performanee  is  shown  as  Figure  26.  A  100  node  run  took 
approximately  368  seeonds.  It  should  be  noted  that  the  optimal  eontrols  for  the  two  burn 
maneuver  oeeur  as  predieted  by  Lawden.  Again  as  diseussed  previously,  differenees 
between  these  results  and  true  impulsive  behavior  ean  be  explained  by  the  diserete  time 
nature  of  this  modeling. 
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Figure  24  State  and  control  histories  for  semimajor  axis  change 


Figure  25  Mass  flow  for  semimajor  axis  change 
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Figure  26  Semimajor  axis  change  transfer  orbit 


3,  Semimajor  Axis  Change  Optimal  Behavior 


a.  Semimajor  Axis  Change  Feasibility 

Feasibility  was  once  again  shown  by  independently  propagating  the  initial 
conditions  through  the  equations  of  motion  using  interpolated  controls.  A  graphical 
representation  of  the  result  is  shown  as  Figure  27.  We  can  conclude  that  the  result  is 
within  the  feasible  set  of  solutions. 
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Figure  27  Plots  of  semimajor  axis  change  model  performance  vs.  ODE45 

propagation 


b.  Semimajor  Axis  Change  Hamiltonian  Behavior 
Graphical  representation  of  the  Hamiltonian  during  the  performance 
period  is  shown  as  Figure  28.  It  can  be  noted  that  the  Hamiltonian  exhibits  a  near 
constant  zero  value  for  the  majority  of  the  run.  By  this  test  condition,  we  can  conclude 
that  optimal  control  is  being  performed. 


Figure  28  Hamiltonian  behavior  for  semimajor  axis  change 
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c.  Semimajor  Axis  Change  Lagrangian  Behavior 
Graphical  representation  of  the  stationary  Lagrangian  during  the 
performance  period  is  shown  as  Figure  29  for  all  six  controls.  We  again  observe  that  the 
minimized  Lagrangian  exhibits  a  near  zero  value  for  the  majority  of  the  run.  By  this  test 
condition,  we  can  conclude  that  optimal  control  is  being  performed. 


Hamiltonian  Minimization 


Figure  29  Semimajor  axis  change  minimized  Lagrangian  behavior 

d.  Semimajor  Axis  Change  KKT  Complementarity  Condition 
Graphical  representation  of  the  switching  structure  for  each  of  the 
control/dual  pairs  is  shown  as  Figure  30.  This  switehing  generally  follows  the  eriteria 
described  by  equation  (4.29),  and  therefore  we  can  conclude  that  optimal  control  is  being 
performed.  It  can  be  observed  that  the  only  eontrols  being  exercised  are  by  the  +  and  - 1 
thrusters,  Tn  and  Tti-  Although  not  shown  again,  expansion  of  these  plots  again 
illustrates  scaling  and  tolerance  issues  noted  previously. 
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Figure  30  Switching  behavior  for  semimajor  axis  change 


4.  Semimajor  Axis  Change  Scenario  Conclusions 

Based  on  the  results  of  this  scenario,  we  can  again  conclude  that  optimal  control 
is  being  performed.  Although  we  have  already  shown  this  for  a  coplanar  transfer  through 
the  Hohmann  scenario,  this  scenario  added  complexities  in  non-zero  eccentric  orbits,  and 
provided  additional  “feel  good”  data  for  the  model  against  a  known  result. 

The  semimajor  axis  change  actually  has  two  known  optimal  solutions  [Chobotov 
(2002)],  one  involving  two  burns  as  has  been  shown,  and  one  involving  a  single 
impulsive  burn  at  either  the  upper  or  lower  intersection  points  of  the  initial  and  terminal 
orbits.  It  should  be  noted  that  given  the  right  conditions,  the  model  can  also  approximate 
the  single  burn  maneuver.  This  solution  is  shown  as  Figure  31,  and  it  meets  all  tests  for 
feasibility  and  optimality.  Again,  it  does  not  exactly  match  the  classical  impulsive 
maneuver,  but  resemblance  to  this  maneuver  is  beginning  to  be  achieved,  and  differences 
can  be  attributed  to  the  particular  constraints  imposed. 
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Figure  31  Single  burn  semimajor  axis  change  maneuver 


E.  MINIMUM  TIME  TRANSFER  SCENARIO 

For  the  final  validation,  we  decided  to  use  the  model  to  perform  optimal  controls 
for  a  minimum  time  problem.  This  involved  some  additional  changes  to  the  model 
parameters,  since  the  other  problems  (and  the  ultimate  aim  of  the  research)  all  involve 
minimizing  fuel.  Flowever,  by  using  the  model  in  a  slightly  different  context,  it  was 
hoped  that  final  validation  of  the  equations  of  motion  could  be  achieved  independent  of 
the  type  of  application. 

This  problem  was  modeled  to  replicate  a  minimum  time  problem  as  reproduced  in 
the  text  by  Bryson  (1999),  but  first  documented  by  Moyer  and  Pinkham  (1964). 

1,  Minimum  Time  Transfer  Input  Parameters 

The  first  major  change  relative  to  previous  validations  was  a  difference  in  the  cost 
functional.  Since  the  goal  of  the  minimum  time  problem  is  obviously  to  execute  a 
maneuver  in  the  minimum  amount  of  time  possible,  fuel  consumption  was  no  longer  to 
be  considered  an  issue.  The  cost  index  was  changed  to  reflect  this,  from  final  mass  to 
final  time. 
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The  input  parameters  for  this  problem  are  based  upon  Bryson’s  textbook 
formulation  (pp.  126-128).  The  unusual  values  of  the  numbers  are  to  simulate  the 
geometry  of  an  Earth  to  Mars  transfer  as  elosely  as  possible,  and  therefore  the  initial 
eonstant  values  for  the  model  as  used  previously  were  modified  to  simulate  Bryson’s 
numerieal  solution.  The  eonstant  fp  was  ehanged  to  492.4567  seeonds  to  mateh  Bryson’s 
eonditions  for  this  problem.  All  other  eonstants  were  left  the  same.  It  should  be  noted, 
sinee  the  problem  as  ealeulated  is  sealed  to  make  the  gravitational  eonstant  equal  to  1 .0,  it 
is  not  neeessary  to  ehange  gravity  to  a  sun  centered  system  for  this  validation. 

The  Keplerian  inputs,  per  Bryson,  were  provided  as  follows: 


a.  =  6378.1  km 

Qj-  =9718.2  km 

e.  =  0.0 

II 

o 

o 

i  =  35  *  rad 

180 

if  =35*-^  rad 
^  180 

Q  =  45  rad 

‘  180 

Q  =  45  *  rad 

^  180 

CO  =  45  *  rad 

'  180 

o.  =  45  *  rad 

^  180 

K  =  0  *  rad 

'  180 

Vf  =235*-^  rad 
^  180 

2.  Minimum  Time  Transfer  Performance 

Given  only  the  above  end  point  conditions,  and  subject  to  the  cost  and  dynamic 
equations,  the  following  performance  was  recorded.  Graphical  representations  of  the 
state  and  control  histories  are  shown  as  Figure  32.  Mass  flow  rate  is  shown  as  Figure  33. 
A  pictorial  plot  of  orbital  performance  is  shown  as  Figure  34.  A  100  node  run  took 
approximately  122  seconds. 
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Figure  32  State  and  control  histories  for  minimum  time  transfer 
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Figure  34  Minimum  time  transfer.  Note  period  of  zero  net  thrust 

An  important  revelation  was  made  through  this  performanee.  It  can  be  readily 
observed  that  this  behavior  is  similar  in  some  respects  to  that  documented  by  Bryson,  but 
it  can  not  be  considered  an  extrapolated  result  as  could  the  other  validation  scenarios. 
Although  the  thrusting  profile  appears  to  spiral  similar  to  what  is  predicted  by  Bryson’s 
work  (Figure  34),  a  significant  portion  of  the  orbit  transfer  undergoes  an  unexpected 
period  of  zero  net  thrust.  This  is  due  a  major  difference  between  Bryson’s  model  and 
ours,  namely  his  model’s  use  of  a  single  gimbaled  thruster  with  constant  thrust  vs.  our  six 
fixed  thruster  model.  Still,  intuition  would  have  us  believe  that  a  constant  thrusting 
profile  would  be  necessary  for  a  minimum  time  maneuver.  This  apparently  troublesome 
point  will  turn  out  to  be  very  significant,  and  will  be  discussed  later  in  the  section. 

3,  Minimum  Time  Transfer  Optimal  Behavior 
a.  Minimum  Time  Transfer  Feasibility 

Feasibility  was  once  again  shown  by  independently  propagating  the  initial 
conditions  through  the  equations  of  motion  using  interpolated  controls.  A  graphical 
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representation  of  the  result  is  shown  as  Figure  35.  From  this,  we  can  conclude  that  the 
result  is  within  the  feasible  set  of  optimal  solutions. 


Figure  35  Plots  of  minimum  time  transfer  model  performance  vs.  ODE45 

propagation 


b.  Minimum  Time  Transfer  Hamiltonian  Behavior 

Graphical  representation  of  the  Hamiltonian  during  the  performance 
period  is  shown  as  Figure  36.  It  can  be  noted  that  the  Hamiltonian  exhibits  a  near 
constant  value  for  the  majority  of  the  run. 
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Figure  36  Hamiltonian  behavior  for  minimum  time  transfer 
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It  is  observed  that  the  constant  value  of  the  Hamiltonian  for  this  case  is  -1 
vs.  the  zero  value  observed  for  the  minimum  fuel  scenarios.  This  can  be  expected  by 
revisiting  the  Hamiltonian  evolution  and  value  condition,  similar  to  the  development 
carried  out  in  equations  (4.19)  through  (4.21).  For  the  minimum  time  case,  the 
Hamiltonian  can  be  written  as: 

//  =  t,+;i,(d)+^^(/;)+4^(P2)+4.(ei)+\(4)+4(^)+^M(^)  (5-7) 

The  Hamiltonian  Evolution  then  indicates  that  a  constant  value  must  be  maintained  for 
optimality,  because: 


dH  _dH 
dt  dt 

The  Hamiltonian  Value  condition  can  then  be  derived  as: 


(5.8) 


BE 

dt^ 


E  =  E  +  K  e. 


■tf+Va 


+  Vpep  +Vpep 


So: 


BE 
Bt  ^ 


-1 


(5.9) 


Combining  the  Hamiltonian  Value  with  the  Hamiltonian  Evolution 
indicates  that  a  value  of  -1  must  be  maintained  by  the  Hamiltonian  during  the  optimal 
control  span.  By  this  test  condition  and  our  observations  then,  we  can  conclude  that 
optimal  control  is  being  performed  during  this  scenario. 


c.  Minimum  Time  Transfer  Lagrangian  Behavior 
Graphical  representation  of  the  stationary  Eagrangian  during  the 
performance  period  is  shown  as  Eigure  37  for  all  six  controls.  We  again  observe  that  the 
stationary  Eagrangian  exhibits  a  near  zero  value  for  the  majority  of  the  run.  This 

particular  scenario  appears  to  exhibit  some  difficulty  maintaining  optimal  control  at 
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periods  in  the  vicinity  of  the  large  spiral  direction  shift  where  thrust  is  momentarily 
directed  in  a  non-optimal  manner,  however  optimality  is  quickly  regained.  By  this  test 
condition,  we  can  conclude  that  optimal  control  is  being  performed  for  the  majority  of  the 
run. 
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Figure  37  Minimum  time  transfer  minimized Lagrangian  behavior 


d.  Minimum  Time  Transfer  KKT  Complementarity  Condition 

Graphical  representation  of  the  switching  structure  for  each  of  the 
control/dual  pairs  is  shown  as  Figure  38.  This  switching  follows  the  criteria  described  by 
equation  (4.29),  and  therefore  we  can  conclude  that  optimal  control  is  being  performed. 
A  detailed  view  of  this  switching  on  a  tighter  scale  for  two  of  the  controls  (Figure  39) 
indicates  that  bang-bang  control  is  clearly  achieved.  This  is  the  type  of  result  sought  but 
not  observed  earlier  in  the  minimum  fuel  scenarios.  We  can  note  in  this  case  that  the 
costate  does  not  remain  near  zero,  but  instead  passes  cleanly  and  quickly  through  the  zero 
crossing,  offering  some  credence  to  the  scaling  and  tolerance  issues  noted  earlier. 
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Figure  38  Switching  structure  for  minimum  time  transfer 


Detailed  Complementarity  Switching  Structure 


Figure  39  Switching  detail  for  minimum  time  transfer,  T,i  and  Tt2 
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4,  Minimum  Time  Transfer  Scenario  Conclusions 

From  the  observations  of  the  preeeding  seetions,  we  ean  again  eonelude  that 
optimal  eontrol  is  being  performed.  However,  as  we  have  noted  previously,  the  thrusting 
profile  of  the  maneuver  does  not  instantly  meet  our  intuitive  expeetations.  Despite  this 
initial  observation,  upon  eloser  inspeetion  and  analysis  we  find  that  the  model  is  aetually 
performing  remarkably. 

We  have  already  drawn  attention  to  the  period  of  zero  net  thrust  observed  in  the 
plot  shown  in  Figure  34.  This  was  unexpeeted  based  on  what  we  were  predisposed  to 
look  for  based  on  Bryson’s  work.  However,  the  differenee  between  six  thrusters  and  one 
(as  modeled  in  Bryson’s  text)  turns  out  to  be  more  surprising  than  we  would  believe.  If 
we  look  eloser  at  the  eontrol  history  shown  in  Figure  32,  we  now  note  two  subtle 
eharaeteristies.  First,  the  positive  and  negative  normal  thrusters  are  firing  at  maximum 
eapaeity,  equal  in  magnitude  and  opposite  in  direetion,  for  the  duration  of  the  seenario. 
The  net  effeet  of  these  firings  is  zero  net  thrust  out  of  the  orbital  plane.  Seeond,  during  a 
short  period  of  time  during  the  seenario  (approximately  TU  1.4  -  1.7),  all  thrusters  are 
firing  at  maximum  eapaeity,  with  eaeh  pair  equal  in  magnitude  and  opposite  in  direetion, 
resulting  in  zero  net  thrust  during  that  short  period. 

The  reason  for  this  behavior  is  explained  simply:  the  model  is  predieting  that 
optimal  eontrol  to  aehieve  a  minimum  time  maneuver  given  the  provided  end  point 
eonditions  (ineluding  the  initial  fuel  load)  is  to  spend  energy  on  zero  net  thrust  maneuvers 
in  order  to  dump  fuel.  By  doing  so,  the  satellite  lightens  its  mass  in  order  to  aehieve  the 
maximum  aeeeleration  possible,  an  obvious  taetie  for  any  maximum  speed  maneuver. 
The  model  ealeulates  how  mueh  fuel  it  needs  to  perform  the  maneuver,  and  offloads  the 
rest  in  order  to  speed  up  exeeution. 

We  ean  ehallenge  this  assertion  by  eommanding  the  normal  thrusters  off  for  the 
duration  of  the  seenario,  leaving  all  other  conditions  the  same  (note:  this  scenario  meets 
all  feasibility  and  optimality  tests  as  provided  previously,  but  are  not  shown  here  for 
brevity).  By  doing  this,  we  can  observe  in  Figure  40  that  the  transfer  takes  more  time  to 
execute  (2.739  TU  vs.  2.365  with  normal  thrusters  on).  This  result  validates  that  fuel 
dumping  provides  for  a  faster  time  solution.  It  is  also  observed  that  the  system  is  firing  in 


66 


plane  throughout  the  profile,  and  does  not  go  through  a  period  of  full  fuel  dump  (Figure 
41),  more  consistent  with  Bryson’s  results.  The  time  of  the  switch  from  positive  radial  to 
negative  radial  still  does  not  quite  match  the  Bryson  profile  (which  occurs  symmetrically 
at  the  transfer  halfway  point),  offering  further  insight  into  the  differences  between  one 
thruster  and  effectively  four  for  this  case.  This  appears  to  be  a  manifestation  of  the 
difference  between  an  f  norm  based  model  and  our  f  norm  based  model  in  associated 
mass  flow  rate  characteristics  [see  Ross  (2004)]. 


Figure  40  State  and  control  history  for  minimum  time  transfer  (normal  thrusters 

disabled) 


5  However,  if  we  inerease  the  starting  mass  of  the  vehicle,  fuel  dumping  again  occurs.  The  lack  of  fuel 
dumping  in  this  case  is  most  likely  due  to  the  use  of  input  values  scaled  from  the  numerical  results  reported 
by  Bryson. 
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Overall,  minimum  time  performance  and  its  unexpected  controlled  fuel  burn  is  the 
final  validation  of  the  equations  of  motion  and  of  the  coding  of  the  model.  The  observed 
behavior  transcends  the  theoretical  mathematics  of  the  model  formulations  and  crosses 
into  predictable  physics  behavior.  For  this  reason,  it  is  believed  this  scenario  provides 
final  proof  that  the  model  is  working  as  it  should,  and  we  are  ready  to  move  on  to  proving 
our  initial  goal;  a  multi-agent  system. 


F.  CONCLUSIONS 

Through  the  scenarios  shown  in  this  section,  it  is  believed  that  the  dynamic 
equations  of  motion  used  and  the  coding  of  the  MATLAB^m  model  have  been  validated. 
The  model  is  ready  to  perform  complex  three  dimensional  maneuvers  which  can  be 
trusted  as  accurate,  and  optimal  behavior  can  be  demonstrated  by  observing  and 
following  the  tenets  of  the  Minimum  Principle.  Complex  control  scenarios  can  be  input 
and  results  observed  following  the  documented  methods,  and  optimal  control  predictions 
can  be  made  within  the  range  of  feasible  solutions  for  any  desired  transfer.  A 
representative  complex  control  result  is  shown  as  Figure  42,  with  changes  made  to  all  six 

Keplerian  elements  between  the  initial  and  final  conditions.  This  scenario  as  run  again 
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meets  all  feasibility  and  optimality  tests,  although  the  details  are  not  shown  for  sake  of 
brevity.  By  using  this  methodology  as  a  baseline,  it  is  hoped  that  multiple  simultaneous 
seenarios  can  be  run  and  predicted  for  optimal  control. 


Figure  42  Complex  orbit  transfer,  all  classical  elements  modified 
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VI.  MULTI-AGENT  CONTROL 


Having  validated  behavior  of  the  orbital  model  for  three  dimensional  motion,  the 
next  step  is  to  modify  it  for  proeessing  multiple  agents.  The  primary  goal  of  this  research 
was  to  prove  that  this  could  be  done  by  effectively  replicating  a  validated  single  satellite 
problem  so  that  multiple  agents  could  be  solved  simultaneously.  The  simultaneous 
distinction  is  important,  since  relationships  can  exist  between  satellites  which  complicate 
matters  tremendously  if  performed  in  a  serial  and  iterative  fashion. 

The  term  “agent”  is  used  generically  to  refer  to  the  body  being  optimized.  For 
this  research,  an  agent  can  be  considered  synonymous  with  a  satellite,  and  the  terms  will 
be  used  interchangeably  within  this  context. 


A.  MODEL  MODIFICATIONS 

In  order  to  modify  the  model  to  process  multiple  simultaneous  satellites,  a  careful 
act  of  replication  had  to  be  carried  out.  Far  from  a  simple  copy  and  paste  exercise,  great 
care  was  maintained  in  instilling  multiple  sets  of  similar  variables  into  the  new  version. 
Certain  variables,  constants,  and  subroutines  were  used  universally  for  all  agents,  while 
others  were  maintained  as  agent  specific.  However,  once  this  process  was  carried  out,  it 
turned  out  to  be  the  most  difficult  part  of  the  leap  from  single  agent  to  multi  agent 
problems. 

A  new  consideration  in  multi-agent  problems  is  establishing  relationships 
between  agents.  This  is  necessary  to  define  certain  conditions  which  the  user  would  like 
to  maintain.  Foremost  is  collision  avoidance:  by  establishing  a  minimum  safe  distance 
between  satellites,  the  model  can  programmed  so  that  two  bodies  cannot  occupy  the  same 
space  at  the  same  time.  Similarly,  a  maximum  distance  can  be  established.  By  using  the 
two  together,  a  relative  operating  window  can  be  created  for  the  purpose  of  establishing 
formation  patterns. 

Another  new  consideration  for  multi-agent  problems  is  establishing  a  new 
performance  index  which  considers  all  satellites.  Again,  for  this  research  minimum  fuel 
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is  the  driving  consideration,  and  the  new  Bolza  cost  used  is  the  simple  combined 
maximized  final  masses  of  the  agents  (that  is,  J  =  -  Mif-  M2f). 

A  third  consideration  for  multi-agent  problems  is  that  care  must  be  taken  not  to 
overconstrain  the  solution.  This  will  be  discussed  in  greater  detail  following  initial 
demonstration  of  the  model,  which  will  help  understanding  of  the  issue  at  hand. 


B,  TWO-AGENT  MODEL 

To  demonstrate  proof  of  the  multi-agent  concept,  the  validated  single-agent  model 
was  replicated  and  coded  to  create  a  two-agent  version.  This  section  documents  the 
performance  of  that  model  following  the  same  development  and  proof  used  in  the 
validation  chapter. 

The  demonstration  of  this  model  will  be  carried  out  using  a  previously  discussed 
validation  case  using  all  three  dimensions.  The  two  satellites  modeled  will  be  assumed  to 
be  launched  on  a  common  launch  vehicle,  then  undergo  equal  and  opposite  orbit  raise 
and  inclination  change.  Because  the  model  must  consider  safe  separation  distance 
between  vehicles,  the  starting  separation  distance  for  the  two  vehicles  will  be  considered 
shortly  after  tipoff  from  the  final  launch  stage. 

Initial  constants  for  the  model  will  also  be  maintained  as  before,  with  the 
exception  that  they  now  need  to  be  replicated  for  each  vehicle.  That  is  to  say: 

^pi=^p2  =220  s 

^01  =^02  =0-00981  km/s^ 

v.i  =v,2  =f.;,(i/2)*go(i/2)  =3.1392  km/s 

^maxl  =  Tmax2  =  1 0  N 

Aoi  =  Ao2  =  398600.4415  km^  / 

Moi  =Mo2  =3000  kg 

It  is  important  to  note  that  generally  speaking  these  values  do  not  all  have  to  be 
equal  for  both  vehicles.  The  gravitational  acceleration  constant  will  always  be  9.81  m/s 
for  all  vehicles  since  that  is  a  propulsion  system  related  constant  that  is  based  on 
Earthbound  testing  and  definition  of  the  thrusters  [Sutton  and  Biblarz  (2001)],  and 

therefore  will  never  vary  regardless  of  the  system  under  consideration  (provided  ,  of 
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course,  that  these  thrusters  are  manufaetured  on  the  Earth).  The  gravitational  parameter  ju 
is  not  required  to  be  equal  for  all  agents;  however,  it  will  in  all  likelihood  be  equal  for 
most  practieal  applieations.  All  other  constants  can  vary  between  vehicles,  and  the  model 
has  eapaeity  for  these  situations^. 

1.  Two-Agent  Input  Parameters 

For  the  two-agent  model,  two  sets  of  initial  and  final  conditions  must  be 
eonsidered.  Following  the  verifieation  eases  in  ehapter  V,  these  parameters  are  as 
follows: 


a.  Agent  1  (Positive  Inclination  Change) 

The  maneuver  by  satellite  1  can  be  eonsidered  much  like  a  combination 
between  a  Hohmann  transfer  and  an  inclination  change.  The  initial  and  final  boundary 
condition  inputs  to  the  model  are  as  follows: 


a.j  =  9567.2  km 

e-i  =  0.0 

L  =35*-^  rad 
180 

Qj,  =10*-^  rad 
“  180 

CO-,  =225*-^  rad 
180 

K.  =0*-^  rad 
180 


=  9567.2  km 

=  0.0 

if,  =55*-^  rad 
180 

Qf.  =10*-^  rad 
180 


=  225*  ^ 

rad 

180 

=  200*  ^ 

rad 

180 

(6.2) 


b.  Agent  2  (Negative  Inclination  Change) 

The  seeond  agent  inputs  are  identieal  to  the  equations  in  (6.2),  with  the 
exeeption  that  the  final  inclination  changes  -20°  rather  than  +20°.  Additionally,  the  initial 


6  Further,  the  model  has  some  limited  capacity  to  deal  with  differences  in  thrust  between  thrusters  on 
the  same  vehicle,  coded  as  an  efficiency  factor  constant  that  can  be  independently  set  for  each  thruster.  For 
this  thesis,  however,  the  efficiency  setting  will  be  left  at  100%  of  the  maximum  available  thrust. 
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true  anomaly  is  slightly  different  that  for  satellite  1  for  the  purpose  starting  at  a  safe 
keepaway  distanee  for  eollision  avoidanee. 


a,i  =9567.20445  km  =1.5  RE 

=  9567.20445  km  =2.0  RE 

o 

o 

II 

=  0.0 

L  =35*-^  rad 

180 

if.  =15*-^  rad 

180 

Q  =  10*-^  rad 
“  180 

Qf.  =  10*-^  rad 

180 

CO-.  =225*-^  rad 
“  180 

(Of.  =  225  *  rad 

180 

K,  =  .1*-^  rad 

180 

Vf,  =  200  *  rad 

”  180 

2,  Two-Agent  Model  Performance 

State  and  eontrol  history  for  satellite  1  (performing  the  positive  inelination 
ehange)  is  shown  in  Figure  43.  History  for  satellite  2  (performing  the  negative 
inelination  ehange)  is  shown  in  Figure  44.  Mass  flow  for  both  satellites  is  shown  in 
Figure  45.  Orbital  performanee  for  both  vehieles  is  shown  in  Figure  46. 

It  ean  be  observed  that  the  two  independent  satellites  elosely  follow  results 
demonstrated  during  the  model  validation  for  a  Hohmann  transfer  eombined  with  an  out 
of  plane  inelination  ehange  maneuver.  The  individual  maneuvers  are  performed  as  we 
would  expeet  from  previous  maneuvers,  with  thruster  firings  oeeurring  at  perigee  and 
apogee  tangent  points,  as  well  as  at  the  nodal  line  for  inelination  ehange. 
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Satellite  1  State  History 
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Figure  43  Satellite  1  state  and  control  histories 


Figure  44  Satellite  2  state  and  control  histories 
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Satellite  1  Mass  Flow 


Figure  45  Satellite  1  and  2  mass  flow 


Figure  46 


Two  satellite  orbit  transfer  with  inclination  change 
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3. 


Two-Agent  Model  Optimal  Behavior 


a.  Agent  1  Feasibility  and  Optimality  Analysis 
Similar  to  the  methods  shown  during  validation,  feasibility  and  optimality 
for  satellite  1  are  shown  below.  Feasibility  is  shown  by  Figure  47.  Flamiltonian  behavior 
is  shown  in  Figure  48.  The  minimized  Flamiltonian  with  respect  to  controls  is  shown  as 
Figure  49.  Switching  behavior  for  satellite  1  is  shown  as  Figure  50  (note  tolerance  issues 
once  again).  Observing  this  data  and  using  the  same  logic  provided  during  validation,  it 
can  be  concluded  that  satellite  1  is  demonstrating  optimal  control. 


Figure  47  Plots  of  satellite  1  performance  vs.  ODE45  propagation 


77 


Hamiltonian 


Figure  48  Multi-agent  Hamiltonian  behavior  (both  satellites) 


Satellite  1  Hamiltonian  Minimization 


Figure  49  Satellite  1  minimized  Lagrangian  behavior 
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b.  Agent  2  Feasibility  and  Optimality  Analysis 
Similar  to  satellite  1,  feasibility  and  optimality  for  satellite  2  are  shown 
below.  Feasibility  is  shown  by  Figure  51.  Hamiltonian  behavior  is  shown  in  Figure  48. 
This  is  the  same  plot  as  for  satellite  1,  since  a  single  Hamiltonian  represents  the  entire 
multi-agent  problem.  The  minimized  Hamiltonian  with  respect  to  controls  is  shown  as 
Figure  52.  Switching  behavior  for  satellite  1  is  shown  as  Figure  53  (again  observing 
scaling  and  tolerance).  Observing  this  data  and  using  the  same  logic  provided  during 
validation,  it  can  be  concluded  that  satellite  2  is  demonstrating  optimal  control. 
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Figure  51  Plots  of  satellite  2  performance  vs.  ODE45  propagation 


Satellite  2  Hamiltonian  Minimization 


Figure  52  Satellite  2  minimized  Lagrangian  behavior 
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4,  Two-Agent  Model  Overconstraint  Issue 

As  performed,  the  two  agent  model  appears  to  provide  aeeurate  prediction  of 
optimal  behavior.  The  individual  satellites  behave  the  same  as  they  do  would  in  a  single 
satellite  model,  and  defined  inter-agent  relationships  appear  to  be  maintained  as  desired. 
This  scenario  helps  fulfill  the  goal  of  our  research,  and  following  the  work  performed  in 
formulating  and  validating  the  problem  seems  rather  anticlimactic.  However,  at  least  one 
point  here  bears  visitation. 

This  scenario  was  largely  symmetric,  and  because  of  this  both  satellites  were  able 
to  finish  in  their  final  required  positions  at  the  same  time.  However,  this  is  not 
necessarily  true  if  “mirror  image”  cases  are  not  requested  for  the  two  satellites.  If 
concurrent  solutions  are  possible,  a  feasible  locally  optimal  solution  may  be  found  that 
meets  the  criteria  required  (i.e.,  an  overconstrained  solution),  but  is  not  the  best  or  most 
fuel-efficient  solution  possible  if  constraints  were  relaxed.  This  is  the  equivalent  to 
providing  a  good  answer  to  a  bad  question. 
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As  an  example,  a  scenario  was  run  requesting  satellite  1  to  perform  an  inclination 
change  only  (maintaining  the  initial  orbital  radius),  and  satellite  2  to  perform  a  Hohmann 
transfer.  All  constraining  end  conditions  were  maintained  similar  to  the  double 
inclination  case.  Phasing  of  the  two  satellites  within  their  respective  final  orbits  was 
requested  to  be  the  same  (200°),  however  since  satellite  2  terminates  in  a  higher  orbit,  it 
requires  a  greater  amount  of  time  to  get  there.  Performance  of  this  maneuvering  is  shown 
as  Figure  54,  and  all  tests  for  feasibility  and  optimality  are  met  (but  not  shown  for 
brevity). 


Figure  54  Multi-agent  maneuver  involving  time  constraint  difficulty.  Note  satellite 

1  orbit  raise  to  loiter. 

It  can  be  observed  that  with  the  discussed  overconstraint  present,  the  optimal 
solution  resulted  in  satellite  1  raising  its  orbit  some  to  in  effect  “loiter”  and  allow  satellite 
2  time  to  reach  its  final  position,  at  which  time  satellite  1  returns  to  its  lower  orbit  to 
reach  its  final  position  at  the  same  time.  This  does  not  match  what  intuitively  makes 
sense  to  us,  whereby  satellite  1  would  remain  at  the  same  altitude  throughout  the 
maneuver  and  change  inclination  at  the  node,  much  like  the  inclination  change  validation 
case  discussed  in  chapter  V.  The  demonstrated  maneuver  appears  comparatively 
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wasteful  of  fuel,  and  is  referred  to  it  as  a  “false  positive”,  however  it  needs  to  be  noted 
that  the  requirement  that  both  satellites  reaeh  their  final  loeations  at  the  same  time  drives 
this  result,  and  that  this  is  in  faet  the  optimal  solution  meeting  these  requirements  (again, 
a  great  answer  to  a  bad  question). 

A  solution  to  this  issue  is  to  eonstrain  only  one  of  the  vehieles  to  a  final  aehieved 
position,  leaving  the  other  to  reaeh  all  orbital  eonditions  required  exeept  for  phasing. 
Instead,  the  seeond  agent  is  allowed  to  phase  itself  in  its  orbit  wherever  it  ean  while 
aehieving  all  other  eonditions,  ineluding  orbit  and  required  operating  distanee  with 
respeet  to  satellite  1.  In  this  way,  a  kind  of  “master/slave”  relationship  between  the 
vehieles  is  implied,  whereby  satellite  1  dietates  to  some  degree  the  final  positioning  of 
the  eonstellation.  But  even  this  is  not  enough  to  avoid  all  problems,  for  some  eare  must 
be  taken  in  this  ease  to  allow  enough  time  for  satellite  2  to  perform  its  neeessary 
maneuvering,  or  else  another  less  fuel  effieient  maneuver  may  oeeur.  For  example,  using 
the  same  model  (sat  1:  inelination  ehange;  sat  2:  Hohmann  transfer)  if  satellite  1  is  again 
allowed  a  phase  ehange  of  200°  to  perform  its  inelination  ehange,  and  satellite  2  is  to 
perform  a  Hohmann  transfer,  it  ean  be  observed  that  enough  time  is  still  not  being 
allowed  for  this  transfer  to  oeeur,  and  therefore  another  solution  to  the  problem  will  be 
found  whieh  is  less  fuel  effieient  (but  meeting  the  eonditions  requested)  (Figure  55). 
This  is  beeause  a  Hohmann  transfer  must  undergo  a  180°  in-plane  phase  shift  in  addition 
to  a  radial  ehange,  whieh  requires  more  time  to  exeeute  than  does  satellite  I’s  maneuver. 
Sinee  final  position  is  tied  to  satellite  1  (the  “master”)  reaehing  200°,  and  in  this  time 
satellite  2  only  shifts  163.4°,  satellite  2  must  expend  additional  fuel  at  both  ends  of  the 
transfer  to  reaeh  its  proper  orbit  within  the  allotted  time.  We  ean  fix  this  diserepaney  by 
either  setting  a  larger  phase  ehange  for  satellite  1  whieh  allows  satellite  2  to  eomplete  a 
180°  maneuver  (Figure  56),  or  by  setting  the  vehiele  with  the  longer  exeeution  time  as 
the  master  vehiele  (Figure  57).  The  net  fuel  eonsumed  is  the  same  for  both  oases. 
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Figure  55 


Figure  56 


Satellite  2  (slave)  not  provided  adequate  time  to  perform  Hohmann 
transfer,  resulting  in  radial  thrusting  to  meet  final  orbit 


Satellite  1  (master)  given  longer  execution  time  to  facilitate  Hohmann 
transfer  by  satellite  2 
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Orbit  Thrust  History 


sat  1  initial  orbit 
-  sat  1  final  orbit 
«  sat  1  transfer  orbit 
sat  2  initial  orbit 
-  •  •sat  2  final  orbit 
a  sat  2  transfer  orbit 
1  thrust  vector 
•sat  2  thrust  vector 


Vfi=200° 

Vf2=unprescribed 


Figure  57  Satellite  1  (master,  performing  Hohmann  transfer)  sets  constellation 
execution  time,  allowing  adequate  time  for  transfer  and  for  satellite  2  inclination 
change.  Note:  satellite  1&2  designations  switched  from  previous  examples 


5.  Two-Agent  Model  Conclusions 

The  two-agent  model  appears  to  predict  optimal  control  behavior  as  expected. 
Model  performance  meets  all  tests  for  feasibility  and  optimality  previously  used  during 
validation,  and  results  of  individual  agents  within  the  multi-agent  construct  appear  to 
behave  the  same  as  similar  single  agent  validations.  Moreover,  by  observing  the 
behavior  of  the  Hamiltonian  for  the  entire  model,  we  can  conclude  that  the  model  as  a 
whole  is  behaving  in  an  optimal  manner.  Some  care  and  consideration  had  to  be  paid  to 
some  slight  behavioral  differences  between  a  single-agent  model  and  a  two-agent  model, 
and  the  two  agent  model  is  more  sensitive  to  how  data  is  input  due  to  time  based 
relationships  between  the  agents.  However,  once  addressed  the  two-agent  model  appears 
to  perform  as  expected. 

The  performance  exhibited  by  this  model  in  a  scenario  with  a  somewhat  expected 
outcome  provides  proof  that  this  model  can  be  run  for  any  number  of  scenarios,  and  using 
the  Minimum  Principle  as  a  guide  optimal  control  can  be  predicted  with  confidence.  A 
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complex  three  dimensional  two-agent  scenario  is  shown  as  Figure  58,  that  meets  all  tests 
for  feasibility  and  optimality  using  the  methods  outlined  throughout  this  thesis. 


Orbit  Thrust  History 


sat  1  initial  orbit 
sat  1  final  orbit 
sat  1  transfer  orbit 
sat  2  initial  orbit 
sat  2  final  orbit 
sat  2  transfer  orbit 
sat  1  thrust  vector 
sat  2  thrust  vector 


Figure  58  Complex  two-agent  scenario 


C.  CONCLUSIONS 

The  proof  of  concept  for  a  two-agent  model  is  a  critical  step  in  proving  full  multi¬ 
agent  optimization  capability,  for  the  replieation  process  for  three  or  more  agents  is 
essentially  the  same  as  the  replieation  proeess  from  one  agent  to  two.  Given  the 
methodology  of  replieation  for  the  two  agent  model,  there  is  no  reason  that  a  greater 
number  of  agents  ean  not  be  similarly  modeled.  At  the  present  time,  it  is  believed  that 
the  only  limitation  on  the  number  of  agents  possible  for  simultaneous  optimization 
involves  processing  power,  which  continues  to  improve  with  time. 
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VII.  CONCLUSIONS  AND  FUTURE  WORK 


A.  CONCLUSIONS 

Based  on  the  researeh  performed  in  support  of  this  thesis,  the  tools  and  methods 
for  multi-agent  eontrol  optimization  have  been  sueeessfully  proven.  By  earefully 
defining  a  set  of  dynamie  equations  for  a  single  agent  model,  and  thoroughly  validating 
these  equations  as  a  foundation  for  the  overall  predietive  model,  it  has  been  eonelusively 
shown  that  optimizing  control  for  multiple  agents  is  as  simple  as  replicating  the 
successful  single  agent  model. 

The  key  to  solving  these  problems  is  adherence  to  Pontryagin’s  Maximum 
Principle,  and  as  has  been  shown  time  and  again  throughout  this  thesis,  fulfilling  the 
necessary  conditions  of  this  principle  provides  convincing  evidence  that  optimal  solutions 
to  a  given  problem  can  be  achieved.  In  truth,  this  research  offers  but  one  small 
application  of  this  theory,  and  its  real  beauty  is  that  it  can  be  extended  to  a  virtually 
limitless  number  of  problems,  the  only  requirement  being  the  ability  to  be  mathematically 
described  by  dynamic  equations. 

The  portability  of  the  process  through  replication  provides  an  easy  avenue  to 
provide  simultaneous  control  of  as  many  agents  as  desired  within  available  processing 
limitations,  and  there  is  little  doubt  that  the  methods  used  here  can  be  used  toward  this 
end.  Further,  as  processing  limitations  continue  to  deplete  through  the  continued 
fulfillment  of  Moore’s  Law,  larger  and  larger  problems  will  be  made  possible,  and 
solutions  will  be  achievable  at  ever  faster  rates.  Just  within  the  research  span  of  this 
thesis,  a  new  release  of  MATLAB™  brought  about  observed  solution  times  that  were  two 
to  three  times  faster  than  the  previous  version  of  the  software.  When  compared  to  similar 
DIDO  related  optimization  projects  of  only  three  to  four  years  ago,  generally  speaking 
solution  times  have  progressed  from  measurement  in  hours  to  measurement  in  minutes. 
This  will  only  continue  to  improve,  enabling  the  tackling  of  more  and  more  complex 
problems  using  these  methods. 


87 


The  model  and  the  tools  eurrently  have  some  limitations  specifie  to  the  studied 
orbital  problem  that  have  not  been  discussed  in  any  detail  thus  far.  They  will  be  captured 
in  the  following  sections. 


B,  MODEL  ISSUES 

1.  Limitation  on  Number  of  Variables 

As  was  touched  upon  in  chapter  V  of  this  thesis,  the  underlying  number  of 
computations  inherent  in  the  model  can  be  envisioned  as  a  directly  linear  function  of  the 
number  of  discrete  points  in  time  (i.e.,  nodes)  chosen  for  any  given  scenario  (equation 
(5.3)).  The  number  of  nodes  chosen  can  be  simply  equated  as  a  measure  of  achievable 
resolution  in  the  problem  solution,  and  this  has  been  illustrated  to  some  degree  in  the 
discussion  on  bang-bang  control. 

All  of  the  scenarios  that  have  been  demonstrated  through  this  thesis  have  involved 
a  small  number  of  orbits,  usually  less  than  one.  Combining  this  with  a  reasonably  large 
number  of  nodes  has  provided  a  fairly  high  resolution  through  the  orbital  path,  which  has 
resulted  in  getting  a  good  picture  of  the  necessary  controls  to  achieve  a  desired  result.  A 
problem  starts  to  rear  its  head,  however,  as  the  number  of  orbits  is  increased  for  a 
maneuver.  This  can  be  illustrated  by  first  envisioning  a  100  node  solution  for  a  maneuver 
that  is  constrained  to  less  than  one  orbit:  here  100  points  can  describe  the  maneuvers 
within  the  one  orbit  span  adequately  for  most  purposes.  Now  suppose  the  maneuver  span 
is  opened  up  to  include  potentially  25  orbits  using  the  same  100  nodes  spread  across  the 
span.  Now  only  four  points  on  average  can  be  considered  per  orbifi.  If  a  significant 
number  of  thrust  points  are  required  (in  a  continuous  sense),  then  much  of  this  detail  will 
be  lost  in  the  result. 

An  obvious  solution  to  this  problem  is  to  continue  to  increase  the  number  of 
nodes  to  as  high  a  resolution  as  possible,  perhaps  only  at  the  expense  of  computation 
time.  However,  the  computation  engine  used  within  the  model  appears  to  suffer  for  high 

7  This  is  not  precisely  true,  since  the  mathematical  workings  of  the  engine  does  not  spread  the  nodes 
linearly  across  the  span.  Instead,  higher  nodal  concentrations  tend  to  cluster  around  the  endpoints  of  the 
span,  which  can  be  easily  seen  in  any  of  the  plots  used  throughout  the  validation  and  multi-agent  sections 
of  the  thesis. 
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numbers  of  nodes  (on  the  order  of  150  and  greater).  Although  this  engine  is  in  theory 
unlimited,  it  involves  eomputations  of  large  matriees  whieh  begin  to  exhibit  greater 
problems  of  ill-eonditioning.  This  presents  a  fundamental  problem  for  higher  numbers  of 
orbits  whieh  require  large  numbers  of  nodes  for  good  resolution,  espeeially  sinee  a  great 
number  of  orbital  problems  are  spread  over  large  timeframes  for  inereased  fuel 
effieieney.  For  these  eases,  the  only  eurrent  workaround  is  to  be  able  to  divide  the 
problem  into  diserete  legs  whieh  ean  be  independently  run.  However,  help  appears  to  be 
on  the  way. 

The  fundamental  drawbaek  in  the  engine  appears  to  originate  with  the  SNOPT 
algorithms  around  whieh  the  DIDO  optimization  tool  is  wrapped.  The  makers  of  SNOPT 
appear  on  the  verge  of  delivering  a  new  version  whieh  is  mueh  more  highly  eapable, 
reportedly  to  the  point  of  managing  variables  in  the  millions  vs.  the  thousands  eurrently. 
If  this  happens,  it  will  render  this  issue  moot,  and  we  ean  again  fall  baek  to  eomplaining 
about  proeessing  speed. 

2.  True  Retrograde  Orbit 

Despite  the  promise  of  zero  singularities  by  using  the  equinoetial  orbital  elements, 
the  model  does  not  work  for  an  absolute  retrograde  orbit  (i.e.-  inelination  of  180°)  as  of 
the  writing  of  this  thesis.  Mathematieally,  the  issue  is  understood  as  has  been  touehed 
upon  in  diseussion  of  the  “retrograde  faetor”  earlier  in  this  thesis,  however  as  a  relatively 
late  addition  to  this  researeh,  full  eoding  issues  throughout  the  many  files,  funetions,  and 
subroutines  of  the  otherwise  operating  model  were  never  fully  resolved  for  this  ease. 
This  problem  begins  to  manifest  itself  greater  as  inelination  approaehes  180°,  and  only 
signifieantly  beyond  even  179°.  For  this  reason,  the  diffieulty  is  eonsidered  a  eoding 
issue  and  left  largely  as  a  footnote  at  this  point,  for  it  is  believed  solvable  with  some 
signifieant  debugging  time,  however  only  useful  for  a  very  limited  number  of  realistie 
oases  and  perhaps  as  an  aoademio  exeroise.  A  great  number  of  other  oases  ranging  from 
0°  to  greater  than  179°  have  been  performed  with  no  diffieulty.  Perhaps  this  is  a  major 
reason  why  the  subjeot  is  oompletely  ignored  in  many  texts,  as  the  vast  majority  of  oases 
are  performed  without  the  added  mathematioal  oomplexity,  no  matter  how  seemingly 
minor. 
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3,  Scaling  and  Tolerance 

As  was  noted  in  several  places  in  this  thesis,  it  is  believed  that  there  exist  some 
issues  of  scaling  and  tolerance  that  are  unresolved  as  of  this  writing.  The  scaling  issues 
appear  to  manifest  in  the  control  duals  in  minimum  fuel  problems,  and  perhaps  this 
coupled  with  behavior  beginning  to  approach  impulsive  behavior  [Lawden  (1963)]  along 
with  internal  zero-tolerance  settings  provides  results  which  are  not  completely  explained 
at  this  point.  It  is  suggested  that  this  be  investigated  further  in  future  work  on  this  topic. 


C.  FUTURE  WORK 

1.  Extension  to  N  Agents 

The  most  obvious  continuation  of  this  work  should  involve  proof  and 
characterization  of  performance  for  greater  numbers  of  satellites.  It  is  believed  at  this 
point  that  this  exercise  should  be  fairly  painless,  however  no  words  can  substitute  for 
actually  performing  the  task  to  show  that  this  is  true.  Some  unanticipated  issues  arose 
with  the  jump  from  one  to  two  agents,  and  some  of  these  issues  (such  as  master/slave 
relationships  and  collision  avoidance  issues)  become  more  complicated  with  the  addition 
of  agents.  The  variable  limitation  issue  discussed  above  will  likely  play  some  factor  in 
the  number  of  agents  achievable  until  worked  out. 

2,  Scaling  Investigation 

There  appear  to  be  issues  which  manifest  in  some  scenarios  run  in  the  model 
where  the  switching  structure  does  not  appear  to  exhibit  clean  bang-bang  behavior  as 
expected,  but  rather  appears  to  peak  at  relatively  low  number  of  nodes  and  “chatter”  at  a 
high  number  of  nodes.  These  scenarios  exhibit  control  costates  which  remain  relatively 
close  to  zero  through  the  duration  of  the  control  chattering,  lending  to  the  belief  that  an 
issue  of  scaling  or  tolerance  may  exist  within  the  model.  Certain  scenarios,  such  as  the 
minimum  time  scenario  run  in  this  research,  do  not  exhibit  this  phenomenon,  while  most 
of  the  minimum  fuel  scenarios  did  to  some  extent.  An  investigation  into  this  behavior 
should  be  continued 
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3.  Increased  Model  Fidelity 

The  model  in  its  eurrent  form  is  ideal  in  the  sense  that  it  does  not  eonsider 
perturbation  or  drag  effects  on  long  term  orbit.  Since  the  thrust  of  this  thesis  was  aimed 
at  injections  of  multiple  satellites  spawned  from  a  single  launch  vehicle,  and  it  is  believed 
that  these  effects  are  very  minor  in  a  relatively  frequent  thrusting  environment,  they  were 
not  concentrated  on  during  the  proof  of  concept  phase  of  this  research.  However,  now 
that  the  concept  has  been  proven,  addition  of  these  effects  is  necessary,  particularly  if  this 
model  is  to  be  extended  to  other  ideas. 

4.  Sensitivity  Analysis 

A  characterization  of  the  model’s  sensitivity  to  various  inputs  vs.  its  optimality 
signatures  should  be  performed  in  order  to  offer  a  better  understanding  of  effects  that  the 
questions  have  on  the  answers.  Such  an  analysis  would  be  useful  in  helping  the  operator 

determine  how  to  best  provide  inputs  which  ensure  that  optimal  solutions  are  provided  in 

the  most  time  efficient  manner. 

5.  Extension  to  Constellation  /  Formation  Control 

The  model  has  been  built  to  be  as  generic  as  possible,  and  throughout  its 
development  an  end  goal  of  realistic  long  term  orbit  simulation  was  envisioned. 
Although  the  concept  of  trajectory  insertion  provided  a  relevant  issue  on  which  to 
concentrate  effort,  ideas  and  problems  of  formation  flying  are  becoming  more  and  more 
widespread  in  the  literature.  It  is  believed  that  these  tools  and  methods  can  be  adapted 
for  use  in  this  regime,  and  therefore  maximum  and  minimum  operating  distances  between 
agents  have  been  accounted  for  in  this  work  with  this  end  in  mind.  It  is  believed  that  this 
area  offers  perhaps  the  most  significant  long  term  extension  of  this  work,  and  the 
prospects  in  this  field  are  exciting.  However,  the  issues  mentioned  in  the  fidelity 
discussion  above  will  very  likely  have  to  be  worked  in  concert  for  a  truly  useful  tool  to  be 
achieved. 
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6,  Proportional  Fuel  Expenditure 

If  formation  modeling  is  attempted,  great  eonsideration  should  be  given  to 
ensuring  that  fuel  burn  vehiele  to  vehiele  is  performed  in  order  to  ensure  proportional  fuel 
eonsumption  oeeurs.  This  will  help  to  ensure  that  all  vehieles  in  the  eonstellation  deplete 
their  fuel  load  at  the  same  time,  and  making  sure  that  the  eonstellation  mission  remains 
intaet  for  as  long  as  possible.  This  model  aeeounts  for  the  possibility  of  different  vehiele 
masses  and  fuel  loads  aeross  all  agents,  so  installation  of  this  type  of  feature  should 
merely  be  applieation  of  additional  path  eonstraints  or  a  similar  means  of  employment. 

7.  Interface  to  Graphical  Visualizer  and  GUI  Input 

Although  never  aceomplished,  an  interfaee  to  some  graphieal  visualization  tool, 
sueh  as  AGFs  Satellite  Toolkit®  would  be  a  very  useful  addition  to  this  tool. 
Additionally,  a  long  term  objeetive  should  be  made  to  make  this  tool  as  user  friendly  and 
transparent  to  the  operator  as  possible.  Some  kind  of  graphieal  user  interfaee  eould  easily 
be  envisioned  whieh  removes  the  user’s  need  to  understand  the  inner  workings  of  this 
model,  instead  eoneerning  him  or  herself  with  only  inputs  and  outputs  and  not  so  much 
on  how  to  feed  the  model. 
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